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INTRODUCTION 


This  final  report  is  for  work  carried  out  under  Grant  No.  AFOSR  88-0139 
during  the  two  year  period  from  April  1,  1988  to  March  31,  1990. 

For  vehicles  that  operate  at  altitudes  above  80  km  such  as  the  Shuttle,  the  Na¬ 
tional  Aerospace  Plane,  and  the  Aeroassist  Flight  Experiment,  calculations  based  on 
the  Navier-Stokes  equations  lead  to  erroneous  results  because  of  the  very  low  densi¬ 
ties.  An  alternative  approach  to  the  use  of  the  Navier-Stokes  equations,  and  one  that 
is  well  suited  to  treating  rarefied  hypersonic  flows,  is  to  employ  the  Direct  Simulation 
Monte  Carlo  (DSMC)  method,  a  statistical  procedure  that  simulates  a  gas  flow  by 
a  large  collection  of  particles.  However,  a  limitation  of  the  DSMC  method  is  that 
it  does  not  allow  efficient  use  of  vector  architectures  that  are  predominate  in  cur¬ 
rent  supercomputers.  Consequently,  the  computation  speed  is  severely  limited,  thus 
restricting  the  problem  size  one  can  handle  to  principally  one-  and  two-dimensional 
flows. 


The  size  of  a  simulation  may  be  measured  in  terms  of  the  number  of  particles 
employed,  and  it  is  clear  that  more  particles  are  needed  for  a  3D  simulation  than  for  a 
ID  or  2D  simulation,  if  the  same  level  of  accuracy  is  to  be  obtained.  Once  the  number 
of  particles  is  set,  the  number  of  cells  into  which  space  is  divided  is  also  set,  because 
the  average  number  of  particles  per  cell  must  be  greater  than  a  certain  minimum  value 
to  reproduce  the  correct  flow  physics.  Practical  experience  shows  this  lower  limit  to 
be  about  15  particles  per  cell,  assuming  modest  statistical  accuracy.  This  number  is 
far  too  small  to  give  useful  statistics  in  simulations  of  nonsteady  flows,  but  for  steady 
flows  one  may  average  the  results  over  time,  effectively  increasing  the  sample  size  per 
cell  to  a  level  that  gives  acceptable  statistics.  To  model  a  3D  problem,  it  would  be 
desirable  to  have  at  least  500, 000  cells  (space  of  128  x  64  x  64  cells)  at  an  average 
density  of  20  particles  per  cell  for  a  total  of  107  particles.  If  the  gas  consists  of  several 
species,  then  a  still  larger  number  would  be  required. 


3asing  an  estimate  on  a  single-species  gas  and  using  103  time  steps  for  the  re¬ 
quired  time  averaging  (20,000  samples  per  cell  on  average)  we  see  that  1010  particle¬ 
time  steps  of  computational  effort  is  required  to  develop  a  solution.  If  the  computation 
is  to  be  carried  out  in  a  period  of  3  hours,  then  a  performance  of  roughly  1  microsecond 
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per  particle  per  time  step  is  required.  Past  experience  has  shown  that  this  perfor¬ 
mance  is  about  one  to  two  orders  of  magnitude  beyond  the  capability  provided  by 
the  DSMC  method  implemented  on  the  Cray  family  of  computers.  Recognizing  that 
the  ratio  in  computation  time  between  a  straightforward  serial  code  and  an  efficiently 
vectorized  code  is  roughly  a  factor  of  15  to  20  on  the  Cray-2,  it  is  clear  that  one 
must  also  consider  other  algorithmic  improvements  beyond  those  directly  related  to 
vectorization  issues,  if  the  necessary  overall  improvement  in  performance  is  to  be 
attained. 

Our  work  has  focused  on  a  reformulation  of  the  DSMC  method  with  the  objec¬ 
tive  of  designing  a  procedure  that  is  optimized  to  the  vector  architecture  found  on 
the  Cray-2.  In  this  regard,  we  have  developed  a  vectorizable  collision-selection  rule 
to  replace  the  time-counter  algorithm  of  the  DSMC  method  so  that  vector  processing 
can  be  more  efficiently  employed.  This  rule,  which  is  used  to  statistically  control  the 
selection  of  both  the  type  and  number  of  colliding  particle  pairs,  has  been  a  key  step 
in  the  development  of  the  method.  Additionally,  our  work  has  focused  on  finding  a 
better  balance  between  algorithmic  complexity  and  the  total  number  of  particles  em¬ 
ployed  in  a  simulation  so  that  the  overall  performance  of  a  particle  simulation  scheme 
could  be  greatly  improved.  Because  of  the  rather  significant  algorithmic  changes  in¬ 
troduced,  it  was  decided  to  first  study  a  single-species,  ideal  gas,  until  full  verification 
of  the  approach  is  achieved,  prior  to  incorporating  real  gas  effects. 


COMPLETED  WORK 

A  major  effort  in  the  period  was  devoted  to  the  writing  of  a  paper  [1],  which 
was  submitted  for  publication  to  the  Physics  of  Fluids  A.  Because  it  outlines  the 
main  objective  of  our  work  and  fully  describes  the  approach  we  are  taking,  a  copy  of 
the  manuscript  is  included  in  Appendix  A.  This  work  presents  a  detailed  discussion 
of  the  development  of  our  theory  fen  a  new  selection  rule  governing  collisions  in  a 
particle  simulation  of  rarefied  hypers.,  c  flow.  The  new  rule  is  particularly  well 
suited  for  use  on  vector-oriented  computers  and  it  was  specifically  developed  with  the 
objective  of  taking  advantage  of  the  much  higher  computation  speed  that  is  available 
when  running  vectorized  codes.  In  the  paper  we  were  able  to  show  that  shock- 
wave  profiles  predicted  by  our  method  agree  exactly  with  the  corresponding  profiles 


2 


1 


predicted  by  the  Direct  Simulation  Monte  Carlo  (DSMC)  method;  and  comparisons 
are  presented  there  for  both  density  and  temperature,  for  hard  sphere  molecules  and 
Maxwell  molecules,  and  for  shock- wave  Mach  numbers  of  3  and  10. 

In  addition  to  the  profile  comparisons,  the  paper  shows  that  the  principal  equa¬ 
tion  on  which  our  theory  is  based  can  be  related  to  the  time-counter  procedure  used 
in  the  DSMC  method.  This  gives  further  proof  that  the  two  methods  are  entirely 
equivalent.  On  the  other  hand,  the  selection  rule  governing  collisions  in  our  theory 
leads  to  an  algorithm  that  can  be  efficiently  implemented  on  a  vector-oriented  com¬ 
puter,  while  the  time-counter  method  leads  to  a  very  inefficient  implementation.  A 
demonstration  of  the  greater  computational  efficiency  that  can  be  realized  on  the 
Cray-2  supercomputer  with  our  method  is  given  by  three  examples  in  the  paper.  The 
three  cases  selected  were  all  carried  out  as  three-dimensional  simulations,  where  the 
number  of  particles  used  ranged  from  106  to  107  and  the  number  of  cells  used  ranged 
from  3  x  104  to  4  x  105. 

A  second  major  effort  was  the  completion  of  a  Ph.  D.  thesis  by  J.  D.  McDon¬ 
ald  [2].  His  thesis  describes  all  of  the  work  that  was  done  to  develop  the  theory  we 
are  using  and  to  create  the  vectorized  code  for  running  on  the  Cray-2.  His  study,  in 
effect,  completes  our  work  on  the  use  of  a  single-species  gas.  Beyond  the  material 
presented  in  Ref.  [1],  his  thesis  gives  a  detailed  account  of  the  extension  of  the  theory 
we  are  using  to  the  case  of  a  multiple  species  gas,  introduces  an  arbitrary  power  law 
for  molecular  interactions,  describes  the  approach  he  has  developed  to  treat  the  ro¬ 
tational  degrees  cf  freedom  for  a  diatomic  molecule,  introduces  his  proposed  method 
for  handling  vibrational  nonequilibrium,  sets  the  stage  for  the  inclusion  of  chemical 
nonequilibrium  in  our  simulations,  and  gives  the  special  programming  steps  taken 
to  implement  vectorized  code  on  the  Cray-2,  as  well  as  the  additional  programming 
steps  taken  to  speed  up  many  operations  in  general. 

For  example,  he  makes  extensive  use  of  lookup  tables,  as  opposed  to  recalculat¬ 
ing  frequently  used  quantities.  A  special  thermalized  reservoir  corresponding  to  the 
freestream  state  is  used  to  supply  particles  for  input  into  the  freestream,  as  opposed 
to  the  use  of  random  numbers  and  extensive  calculations  to  generate  this  state.  This 
same  reservoir  of  thermalized  particles  is  used  as  a  source  of  samples  for  determining 
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the  division  of  energy  between  vibration,  rotation,  and  translation  for  colliding  par¬ 
ticles,  as  opposed  to  employing  the  Borgnakke-Larsen  method  with  its  heavy  use  of 
exponential  functions.  In  addition,  rather  complete  information  on  the  programming 
used  is  presented  along  with  a  number  of  examples  of  the  computed  results  for  differ¬ 
ent  simulated  conditions.  His  thesis  represents  the  basis  on  which  all  of  our  further 
investigations  will  be  carried  out.  Also,  his  work  has  established  the  approach  we  will 
be  using  to  carry  out  the  vectorization  of  the  work  remaining  to  be  done. 

Different  aspects  of  the  above  work  were  also  reported  in  two  AIAA  papers 
which  were  given  at  the  24th  Thermophysics  Conference  at  Buffalo  in  June,  1989.  A 
paper  by  Woronowicz  and  McDonald  [3]  compares  the  simulated  flow  past  a  wedge 
with  results  from  a  corresponding  experiment.  Although  the  maximum  Reynolds 
number  that  could  be  achieved  in  the  simulation  (for  a  reasonable  computation  time) 
was  below  the  value  reported  for  the  experiment  (3,560  versus  13,500),  density,  tem¬ 
perature,  and  pressure  distributions  were  shown  to  compare  very  well,  while  the 
velocity  profiles  in  the  boundary  layer  were  found  to  be  similar  but  thicker  due  to  the 
smaller  Reynolds  number. 

The  paper  by  Feiereisen  and  McDonald  [4]  is  a  product  of  the  close  collabora¬ 
tion  between  the  work  being  done  by  our  Stanford  group  and  work  by  NASA-Ames 
scientists  in  the  Aerothermodynamics  Branch.  In  this  case,  Feiereisen  (NASA-Ames) 
developed  the  programming  for  creating  an  arbitrary  three-dimensional  body  on  an 
IRIS  Workstation,  which  then  outputs  to  a  file  the  resulting  body  geometry  and 
associated  constants  needed  for  the  boundary  conditions.  This  file  is  then  used  in 
conjunction  with  McDonald’s  particle  simulation  program  to  obtain  a  solution  for 
the  rarefied  hypersonic  flow  about  the  body.  This  collaborative  effort  led  to  the  first 
successful  attempt  to  study  the  full  three-dimensional  flow  about  the  body  shape 
representing  the  Aeroassist  Flight  Experiment  (AFE). 

Shown  in  Figs.  1  and  2  are  the  resulting  solutions  for  the  pressure  and  temper¬ 
ature  distributions  in  the  plane  of  symmetry  and  in  a  downstream  wake  cross-section 
of  the  AFE.  The  blunt  body  shape  was  simulated  accurately  by  the  front  surface, 
but  it  was  terminated  by  a  flat  plane  on  the  aft  side.  The  stair-stepped  outline  seen 
is  a  result  of  the  use  of  a  graphical  polygon-fill  routine  to  create  the  front  surface, 
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but  the  terminating  line  was  not  defined.  These  runs  were  performed  with  9.5  x  106 
particles  in  a  space  of  120  x  60  x  60  cells  for  a  total  of  4  x  105  mesh  cells.  The 
body  diameter  was  44  cell  lengths,  the  upstream  mean- free  path  length  was  1.0  cells, 
and  the  Mach  number  was  35.  Based  on  the  body  diameter,  the  Reynolds  number 
was  2,300.  Good  statistics  were  obtained  with  time  averaging  over  800  time  steps 
using  a  total  of  4.5  hrs.  of  Cray-  2  single  processor  CPU  time.  Because  of  the  large 
data  sets  required,  these  simulations  made  use  of  about  125  MW  of  Cray-2  central 
memory. 


MULTIPLE  SPECIES 

All  of  the  work  reported  above  was  carried  out  for  the  rase  of  a  single-species 
gas.  The  stage  was  reached  recently  where  we  were  confident  that  our  approach  had 
been  adequately  tested  and  that  it  is  basically  correct.  Following  this,  considerable 
work  was  devoted  to  the  creation  of  new  code  for  handling  multiple  species.  This  code 
was  written  in  the  C  programming  language  by  Jeffrey  McDonald,  and  time-critical 
elements  of  the  code  were  being  considered  for  translation  into  assembly  language  for 
the  Cray-2,  just  as  was  done  for  the  single-species  code  that  we  have  been  using.  The 
C  version  of  this  code  was  created  and  tested  before  McDonald  finished  his  thesis  and 
he  was  able  to  give  a  full  discussion  of  his  efforts  there.  Some  of  his  tests  involved  the 
computation  of  the  shock-wave  profile  for  a  mixture  of  two  hard-sphere,  monatomic 
gases;  and  one  of  his  test  results  is  shown  in  Fig.  3,  where  a  comparison  is  made 
between  his  predictions  and  results  obtained  using  the  DSMC  method.  As  can  be 
seen,  the  two  agree  very  nicely  and  McDonald  was  able  to  conclude  that  his  vectorized 
approach  was  working  properly.  Although  McDonald  finished  his  work  as  a  student, 
he  still  participates  in  our  research  efforts  as  he  now  has  a  position  at  NAS  A- Ames, 
where  he  is  associated  with  the  Numerical  Aeronautical  Simulation  (NAS)  Program, 
and  has  been  given  responsibility  for  work  that  is  closely  connected  to  his  previous 
work  as  a  student. 

At  the  start  of  calendar  year  1990,  several  changes  were  made  at  the  NAS 
computation  facility  at  NASA-Ames,  which  we  use  to  carry  out  our  simulations.  The 
change  that  has  become  most  exciting  to  us  was  an  upgrade  in  the  version  of  the  C 
compiler  available  for  use  on  the  Cray- 2  and  the  Cray  YMP.  Once  McDonald  and 
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Fallavollita  started  testing  the  new  compiler  and  set  their  sights  on  continuing  the 
development  of  McDonald’s  code  for  multiple  species,  which  Fallavollita  is  taking  on 
as  his  thesis  topic,  it  became  evident  that  the  new  C  compiler  contains  new  features, 
and  a  new  capability,  which  greatly  reduces  the  amount  of  work  we  will  have  to  do. 

The  new  features  now  provide  the  capability  in  the  C  programming  language 
that  earlier  could  only  be  accessed  after  McDonald  created  considerable  assembly 
language  code,  unique  to  the  Cray-2.  In  fact,  essentially  all  o ''  our  previous  assembly 
language  code  can  now  be  replaced  by  direct  programming  in  C,  assuming  appropriate 
attention  is  given  to  the  use  of  proper  code  constructs.  McDonald’s  tests  have  shown 
that  the  new  compiler  has  all  of  the  capabilities  that  we  need  and  that  it  produces 
code  that  executes  just  as  fast  as  his  earlier  assembly  language  code.  This  means  that 
we  can  write  the  programming  for  multiple  species,  and  for  the  inclusion  of  chemical 
nonequilibrium,  in  a  much  shorter  time  than  would  have  been  possible  with  the  earlier 
version  of  the  C  compiler.  In  addition,  this  code  can  now  be  run  on  either  the  Cray-2 
or  the  Cray  YMP.  Although  the  NAS  Cray  YMP  has  128  MW  of  central  memory 
versus  256  MW  for  the  NAS  Cray-2,  it  is  twice  as  fast  and  represents  a  capability 
which  we  are  interested  in  using.  This  development  represents  a  very  significant  step 
forward  in  the  rate  at  which  we  will  be  able  to  reach  our  goal  of  creating  vectorized 
code  that  includes  chemical  nonequilibrium.  The  first  report  that  will  result  from  this 
new  capability  is  an  AIAA  paper  [5]  to  be  given  by  Feiereisen  and  McDonald  at  the 
5th  Joint  Thermophysics  and  Heat  Transfer  Conference  in  Seattle  in  June,  1990.  This 
paper  will  discuss  our  progress  made  in  dealing  with  multiple  species  and,  hopefully, 
initial  steps  toward  the  inclusion  of  nonequilibrium  effects. 


CHEMICAL  NONEQUILIBRIUM 

Concurrent  with  the  development  of  the  code  for  multiple  species,  Brian  Haas 
has  been  researching  the  problem  of  extending  our  simulation  method  to  include 
the  treatment  of  chemically  reacting  flows,  while  retaining  the  advantages  inherent  in 
vectorization  present  in  the  current  code.  His  primary  concern  has  been  to  develop  the 
theory  so  that  good  compatibility  would  be  found  between  the  needs  of  the  algebra  for 
modeling  chemical  relaxation  and  the  unique  needs  of  vectorization.  In  this  regard  he 
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and  McDonald  have  worked  very  closely  so  that  the  desired  features  would  be  present 
in  the  method  selected  as  a  result  of  his  study. 


One  area  in  which  Haas  has  made  very  nice  progress  is  in  his  treatment  of 
rotational  and  vibrational  nonequilibrium.  Following  McDonald’s  lead,  he  makes 
use  of  the  known  relation  between  continuous  and  a  quantized  two-degree-of-freedom 
systems  to  cause  the  rotational  (continuous)  and  vibrational  (quantized  harmonic 
oscillator)  states  to  reach  a  common  temperature  (equilibrium),  by  mixing  their  en¬ 
ergies  in  a  way  consistent  with  the  fundamental  physical  assumption  of  statistical 
mechanics.  This  assumption  holds  that  two  states  having  the  same  energy  with  no 
degeneracies  are  equally  probable.  The  primary  purpose  in  using  this  approach  is 
that  the  mixing  of  energies  can  be  done  using  mathematics  that  does  not  involve 
extensive  use  of  exponentials  and  transcendental  functions,  as  is  the  case  with  the 
Borgnakke- Larsen  method  —  the  approach  most  frequently  used  in  Bird’s  Direct  Sim¬ 
ulation  Monte  Carlo  method.  This  is  very  important  in  minimizing  the  computation 
time  that  is  taken  in  carrying  out  a  simulation. 

Some  of  the  basic  problems  he  has  addressed  have  related  to  the  question  of 
how  to  treat  recombination,  which  has  not  been  adequately  treated,  and  even  ig¬ 
nored,  in  work  by  others.  They  have  assumed  that  recombination  would  play  a  minor 
role  in  hypersonic  blunt-body  flows  and  dissociation  would  be  the  dominate  feature. 
The  principal  reason  for  this  assumption  is  that  the  theory  for  recombination,  at  the 
molecular  level,  has  not  been  fully  worked  out.  The  difficulty  one  encounters  is  that 
recombination  is  basically  a  three-body  process  and  a  three-body  collision  is  very  dif¬ 
ficult  to  handle  in  a  particle  method,  where  extremely  rapid  computation  is  essential. 
Also,  there  is  no  clear  physics  associated  with  the  way  one  should  distribute  the  total 
energy  of  the  system  under  recombination  in  a  way  that  meets  all  the  requirements 
associated  with  matching  the  temperature  dependence  of  the  recombination  rate  coef¬ 
ficient,  reproducing  the  correct  equilibrium  distribution  for  the  different  species,  and 
establishing  a  common  temperature  among  all  the  species  of  the  mixture. 

In  Haas’  work  he  has  made  careful  use  of  detailed  balance,  which  is  a  very  impor¬ 
tant  theoretical  concept  in  both  equilibrium  and  nonequilibrium  thermodynamics.  In 
addition,  he  has  made  good  use  of  a  very  reasonable  assumption  that  the  total  energy 
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of  the  system  in  recombination  can  be  distributed  among  the  component  energies  by 
increasing  each  of  them  in  proportion  to  the  energies  they  had  before  recombination. 
This  method  is  exact  for  an  equilibrium  gas,  but  is  not  fully  correct  for  a  gas  that  is 
wholly  out  of  equilibrium.  However,  because  many  other  thermalizing  processes  are 
taking  place  at  a  much  faster  rate  than  recombination,  and  the  fact  that  the  fraction 
that  undergoes  recombination  is  very  small,  the  resulting  approximation  is  a  very 
good  one,  and  he  has  shown  that  it  works  extremely  well. 

Besides  the  presentation  in  his  thesis  which  is  expected  to  be  finished  in  calendar 
year  1990,  he  is  in  the  process  of  writing  several  papers  on  his  work.  The  first  of 
these  will  be  an  AIAA  paper  [6]  to  be  given  at  the  5th  Joint  Thermophysics  and  Heat 
Transfer  Conference  in  Seattle  in  June,  1990.  He  has  been  working  diligently  on  it 
and  the  most  recent  draft  of  his  paper  is  presented  in  Appendix  B.  Here,  we  can  refer 
to  his  figures  4  through  8  to  see  how  well  his  model  works  and  to  assess  the  state 
of  his  accomplishments.  Figure  8  shows  that  he  is  able  to  fully  handle  a  5-species 
representation  of  air  at  temperatures  that  correspond  to  re-entry  conditions. 


FLAT  PLATE  BOUNDARY  LAYER 

Michael  Woronowicz  has  been  studying  the  drag  and  heat  transfer  for  a  flat 
plate  boundary  layer  by  researching  the  literature  and  by  carrying  out  simulations 
with  our  single-species  code.  The  purpose  of  his  work  is  to  provide  our  group  with 
the  knowledge  and  experience  needed  so  that  we  will  be  able  to  treat  and  handle 
the  more  complex  boundary  layers  encountered  in  three-dimensional  flows  about  the 
AFE  body.  It  is  our  view  that  we  must  show  that  we  axe  able  to  predict  the  simple 
case  of  a  flat-plate  boundary  layer  quite  well,  before  we  attack  the  more  complex  case 
of  interest  in  our  work. 

He  has  taken  two  very  important  steps  which  place  our  group  in  a  strong  posi¬ 
tion  to  continue  along  our  present  path.  First,  he  has  devised  a  new  way  to  correlate 
all  of  the  existing  experimental  data  so  that  the  data  collapse  onto  a  single  curve, 
instead  of  appeaxing  as  a  family  of  curves  with  a  free  parameter;  and  second,  he  has 
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conducted  numerous  simulations  with  our  code  to  compare  results  with  these  exper¬ 
imental  data,  which  show  that  we  are  able  to  make  predictions  that  compare  with 
experimental  results  very  well. 

Figure  4  is  taken  from  a  paper  that  Woronowicz  has  been  writing.  It  compares 
all  of  the  available  experimental  data  on  flat-plate  drag  for  the  case  of  an  adiabatic 
plate  and  for  three  different  correlation  methods:  a)  use  of  the  viscous  interaction 
parameter;  b)  the  so-called  slip  parameter;  and  c)  the  new  parameter  introduced  by 
Woronowic2  It  is  very  clear  that  the  new  parameter  is  quite  superior  to  the  other  two 
in  correlating  the  experimental  data.  Woronowicz  is  able  to  show  that  this  parameter 
works  equally  well  for  the  case  of  a  highly  cooled  plate,  as  well  as  for  the  case  of  heat 
transfer  to  a  cold  flat  plate. 

In  Figs.  5  and  6  we  show  comparisons  between  his  simulated  results  and  the 
results  of  experiment,  for  the  case  of  a  cold  flat  plate  and  the  case  of  an  adiabatic,  or 
hot,  flat  plate.  His  new  correlation  parameter  is  used  in  these  two  plots  and  it  is  clear 
that  it  is  working  extremely  well.  At  present,  more  runs  have  been  conducted  for  the 
case  of  the  cold  wall,  and  therefore,  Fig.  5  appears  more  complete  than  Fig.  6.  The 
figures  show  that  the  simulations  are  giving  very  good  results  over  a  wide  range  of 
conditions  and  that  they  blend  in  nicely  with  the  data  in  the  near-continuum  regime. 

Because  of  the  large  number  of  particles  that  can  be  used  in  our  simulations, 
we  are  able  to  obtain  nice  resolution  and  striking  pictures  of  the  simulated  flows.  A 
typical  example  is  shown  in  Fig.  7  where  the  Mach  number  distribution  in  the  flow 
near  the  flat  plate  is  displayed.  It  is  of  great  interest  to  observe  the  velocity  slip  that 
occurs  near  the  leading  edge  of  the  plate.  This  is  strictly  a  macroscopic  phenomena 
that  occurs  in  a  rarefied  flow,  as  the  boundary  condition  used  for  the  individual 
reflected  particles  gives  these  particles  random  thermal  components  of  velocity,  at 
the  plate  temperature,  but  no  mean  motion. 


EFFECT  OF  SCALE 

An  important  effect  that  has  not  been  studied  is  the  connection  between  com¬ 
putation  time  and  the  size  of  a  simulation.  One  of  the  most  useful  properties  of  a 
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particle  code  which  is  based  on  statistical  sampling  is  that  the  computational  cost  is 
proportional  to  the  number  of  particles  and  the  number  of  steps;  that  is,  the  same 
computational  cost  is  incurred  when  doubling  the  number  of  particles  as  when  the 
number  of  time  steps  is  doubled.  However,  if  doubling  the  number  of  particles  leads 
to  a  better  physical  representation  of  a  gas  flow  than  doubling  the  duration  of  the 
averaging,  then  it  is  clear  which  arrangement  one  would  choose,  provided  sufficient 
computer  memory  is  available. 

If  one  assumes  that  the  theory  for  statistically  independent  random  variables 
applies  to  a  particle  simulation,  then  one  concludes  that  the  fractional  error,  e,  asso¬ 
ciated  with  a  macroscopic  variable,  depends  only  on  the  size  of  the  statistical  sample. 
Therefore,  it  is  given  by  e  =  l/y/NT,  where  N  is  the  number  of  particles  in  a  cell 
and  T  is  the  number  of  times  steps  used  in  time  averaging.  Because  computation  cost 
depends  directly  on  the  product  NT ,  one  concludes  that  the  fractional  error  remains 
fixed,  if  cost  is  held  fixed,  independent  of  the  size  of  the  simulation,  N.  On  this  basis, 
it  is  obvious  that  one  would  choose  to  use  a  longer  time  average  because  it  would  tie 
up  less  of  the  computer  resources.  On  the  other  hand,  the  situation  may  be  quite 
different  if  statistical  independence  does  not  fully  apply. 

Figure  8  shows  the  results  of  a  very  preliminary  test  using  our  single-species 
code  developed  for  the  Cray-2.  The  test  problem  studied  was  the  blunt-body  flow 
produced  by  a  two-dimensional  flat  plate  in  a  Mach  8  flow  of  an  ideal,  diatomic, 
hard-sphere  gas,  where  the  Knudsen  number,  based  on  the  plate  width,  was  0.1  and 
the  plate  temperature  was  held  constant  at  the  freestream  value.  The  error  was 
computed  as  the  mean  rms  error  of  the  pressure  in  a  cell  divided  by  the  stagnation 
point  pressure.  Fortunately,  the  absolute  error  can  be  determined  simply  from  a  series 
of  comparisons  of  the  relative  errors  between  the  different  runs,  without  having  to 
know  the  exact  solution.  The  curve  shown  is  for  a  constant  value  of  NT;  and  we 
see  that  the  fractional  error  decreases  dramatically  as  the  number  of  particles  in  the 
simulation  is  increased,  even  though  the  total  computation  cost  is  held  fixed.  The 
figure  indicates  that  the  nonlinear  physics  present  in  the  simulation  requires  that,  in 
order  to  reach  a  level  of  error  of  0.1  percent,  it  is  better  to  employ  a  simulation  having 
an  average  density  of  100  particles  per  cell  than  to  use  10  particles  per  cell  together 
with  a  longer  time  average.  In  fact,  the  preliminary  data  indicate  that,  for  a  given 
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level  of  error,  a  certain  minimum  number  density  is  needed,  regardless  of  the  duration 
of  the  time  average  used.  This  clearly  helps  to  explain  how  to  minimize  computation 
cost  in  a  particle  simulation.  This  work  has  been  carried  out  by  Michael  Fallovollita 
and  the  results  of  a  more  complete  study  are  to  be  published  soon. 


CONTINUING  WORK 

Over  the  past  four  years,  we  have  enjoyed  a  very  close  working  relationship  with 
the  Aerothermodynamics  Branch  at  NAS  A- Ames  Research  Center,  in  particular  with 
Dr.  G.  S.  Deiwert  (Branch  Chief)  and  Dr.  W.  J.  Feiereisen  (Assistant  Branch  Chief), 
and  have  received  support  from  their  branch  in  a  number  of  ways:  grant  funds  to  sup¬ 
port  students,  use  of  their  workstations  and  facilities  by  our  students,  collaborative 
effort  between  individuals  in  the  Stanford  and  Ames  groups,  and  access  to  the  com¬ 
putation  facilities  available  through  the  Numerical  Aerodynamic  Simulation  (NAS) 
Program  at  Ames.  At  the  present  time  our  Stanford  group  consists  of  four  Ph.D. 
students,  and  three  scientists  from  NASA-Ames  regularly  attend  our  weekly  research 
meetings  (Feiereisen,  McDonald  and  Boyd). 

The  collaborative  effort  between  our  respective  groups  was  specifically  planned 
to  consist  of  an  arrangement  where  our  students  would  focus  principally  on  theoretical 
questions  and  on  computations  that  exercise  and  test  new  ideas,  while  the  Ames  effort 
would  give  greater  focus  to  the  solution  of  practical  problems,  employing  the  new 
computational  capability  developed.  The  main  thrust  of  our  continuing  research  is 
to  push  forward  on  the  following  four  fronts: 

1)  to  include  chemical  kinetics  in  our  simulations  so  that  thermochemical  nonequi¬ 
librium  among  several  different  species  may  be  handled; 

2)  to  further  develop  our  graphics  display  capability,  along  with  suitable  data 
analysis,  to  handle  the  huge  data  sets  we  encounter  in  our  simulations; 

3)  to  study  rarefied  hypersonic  boundary  layers  on  flat  plates  to  guide  our  under¬ 
standing  of  boundary  layers  on  complex  three-dimensional  bodies;  and 
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4)  to  apply  these  capabilities  to  the  solution  of  practical  problems  of  rarefied,  real- 
gas,  hypersonic,  three-dimensional  flows  about  aerodynamic  bodies  such  as  the 
AFE. 

As  the  different  efforts  described  above  converge  to  a  single  working  code,  we  will 
turn  our  attention  to  the  solution  of  several  practical  problems  involving  real  vehicles 
and  real-gas  flows.  Here,  our  close  collaboration  with  the  NASA-Ames  Aerothermo- 
dynamics  group  will  be  extremely  useful,  as  the  work  will  require  the  specification  of 
reaction  rates  for  different  species,  the  use  of  realistic  body  geometries,  questions  of 
physical  processes  associated  with  gas-surface  interactions,  conditions  of  heat  transfer 
at  the  body  surface,  and  the  altitude-Mach  number  flight  envelope.  One  of  the  body 
geometries  we  plan  to  investigate  thoroughly  is  the  AFE  body.  These  will  consist  of 
fully  three-dimensional  simulations  employing  over  107  particles,  0.5  x  106  cells,  and 
multiple  species.  The  anticipated  results  ought  to  be  very  exciting  and  we  are  looking 
forward  to  the  point  in  our  work  where  we  will  be  in  a  position  to  study  data  from 
such  runs. 
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Fig.  1 .  Pressure  distribution  in  the  central  plane  of  the  AFE  body  and  in  a  downstream 
wake  cross-section  of  the  flow.  The  simulation  models  the  AFE  body  at  an 
altitude  of  90  km  and  a  Mach  number  of  35.  The  gas  is  nitrogen,  treated  as 
ideal. 
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Fig.  2.  Temperature  distribution  in  the  central  plane  of  the  AFE  body  and  in  a  down¬ 
stream  wake  cross-section  of  the  flow.  The  simulation  models  the  AFE  body  at 
an  altitude  of  90  km  and  a  Mach  number  of  35.  The  gas  is  nitrogen,  treated  as 
ideal. 
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Normalized  Density 


Fig.  3.  Density  distribution  in  a  shock  wave  for  a  hard-sphere,  monatomic,  binary 
mixture  of  gases  and  a  Mach  number  of  10.  The  heavy  species  is  15  times  as 
massive  as  the  light  species  and  is  present  in  a  concentration  of  0.20  of  that  of 
the  light  species. 
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Fig.  4. 


Drag  coefficient  for  a  flat  plate  for  three  different  hypersonic  correlation  parame¬ 
ters:  a)  use  of  the  viscous  interaction  parameter;  b)  the  so-called  slip  parameter; 
and  c)  the  new  parameter  introduced  by  Woronowicz. 
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Fig.  5.  Comparison  of  simulated  results  versus  experiment  for  the  drag  coefficient  on  a 
cold  flat  plate.  Simulations  designated  by  the  solid  symbol. 
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Simulation  vs.  Experiment,  hot  wall  data 


Fig.  6.  Comparison  of  simulated  results  versus  experiment  for  the  drag  coefficient  on 
an  adiabatic  flat  plate.  Simulations  designated  by  the  solid  symbol. 
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Fig.  7.  Mach  number  distribution  in  the  flow  near  a  flat  plate.  Simulations 
carried  out  for  a  freestream  Mach  number  of  21.8,  plate  temperature 
equal  to  the  freestream  temperature,  and  a  hard-sphere,  ideal,  diatomic 
gas. 
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Fig.  8.  Fractional  error,  e,  for  the  pressure  variable,  versus  the  average  number  of 
particles  per  cell,  N,  for  a  fixed  compuational  effort.  Curve  shown  is  for  the 
computational  cost  NT  =  2xl05. 
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APPENDIX  A 


A  collision-selection  rule  for  a  particle  simulation  method 
suited  to  vector  computers 

D.  Baganoff  and  J.  D.  McDonald 
Department  of  Aeronautics  and  Astronautics 
Stanford  University,  Stanford,  CA  94305 


Abstract 

A  theory  is  developed  for  a  selection  rule  governing  collisions  in  a  particle  simulation 
of  rarefied  gas-dynamic  flows.  The  selection  rule  leads  to  an  algorithmic  form  that  is  highly 
compatible  with  fine  grain  parallel  decomposition,  allowing  for  very  efficient  utilization  of 
supercomputers  having  vector  or  massively  parallel  SIMD  architectures  (single  instruction 
multiple  data).  A  comparison  of  shock- wave  profiles  obtained  using  both  the  selection  rule 
and  Bird’s  Direct  Simulation  Monte  Carlo  (DSMC)  method  show  excellent  agreement.  This 
serves  to  establish  the  validity  of  the  method,  as  the  DSMC  method  is  known  to  compare 
well  with  experimentally  determined  shock-wave  profiles.  In  addition,  the  equation  on  which 
the  selection  rule  is  based  is  shown  to  be  directly  related  to  the  time-counter  procedure  in 
the  DSMC  method,  further  establishing  their  equivalence.  The  results  of  several  example 
simulations  of  representative  rarefied  flows  are  presented,  for  which  the  number  of  particles 
used  ranged  from  106  to  107,  demonstrating  the  greatly  improved  computational  efficiency 
of  the  method. 
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I.  INTRODUCTION 


The  most  widely  used  particle  method  for  simulating  a  rarefied  gas-dynamic  flow  is  the 
Direct  Simulation  Monte  Carlo  (DSMC)  method.1,2  This  method  was  introduced  by  G.  A. 
Bird  in  the  1960s  and  has  been  advanced  to  the  state  where  it  is  reliable,  quite  accurate  and 
can  handle  rather  general  conditions  of  gas  complexity  and  flow  geometry.  The  main  problem 
that  arises  in  applying  the  method  is  that  the  principal  algorithm  employed  contains  logic 
which  leads  to  conditional  branching  of  a  type  that  does  not  allow  one  to  make  effective  use 
of  supercomputers  having  a  vector  architecture,  such  as  the  Cray-2,  or  a  massively  parallel 
SIMD  architecture  (single  instruction  multiple  data)  such  as  the  Connection  Machine. 

Because  the  execution  speed  of  an  efficiently  vectorized  code  is  over  an  order  of  mag¬ 
nitude  faster  on  a  vector  oriented  machine  such  as  the  Cray-2  than  the  speed  of  the  same 
code  when  limited  to  scalar  execution,  it  is  important  to  explore  methods  that  fully  utilize 
these  features.  In  addition,  many  of  the  projections  estimating  future  increase  in  compu¬ 
tational  capability  are  principally  based  on  extensive  use  of  parallelism,  requiring  scientific 
applications  to  be  properly  structured  to  use  such  architectures  effectively. 

It  is  well  known  that  Monte  Carlo  schemes  in  general  do  not  vectorize  well,  as  this 
problem  has  been  addressed  a  number  of  times,  particularly  in  the  area  of  plasma  physics. 
In  a  recent  fluid  mechanical  study,  Ploss  vectorized  a  modified  Nanbu  algorithm,  a  simulation 
method  for  solving  the  Boltzmann  equation,  and  concluded  that  the  predominantly  scalar 
Bird  scheme  was  still  the  faster  approach.3  However,  in  light  of  the  fact  that  over  an  order 
of  magnitude  decrease  in  computation  time  could  be  realized  if  the  DSMC  approach  could 
be  restructured  to  operate  in  a  vector  fashion,  there  is  a  strong  motivation  to  seek  ways 
to  reformulate  the  algorithms  utilized  in  the  DSMC  method  which  are  presently  limited  to 
scalar  execution. 

To  identify  the  elements  that  slow  computation  in  the  DSMC  method  and  to  under¬ 
stand  the  kind  of  changes  required,  we  first  outline  the  general  steps  that  characterize  a 
particle  method  and  then  focus  the  discussion  on  the  particular  step  which  represents  Bird’s 
time-counter  procedure  and  the  problems  associated  with  its  use.  This  is  followed  by  the 
development  of  a  theory,  leading  to  an  algorithmic  structure  for  a  selection  rule  governing 
collisions  that  is  highly  compatible  with  fine  grain  parallel  decomposition.  A  comparison  of 
shock-wave  profiles  obtained  using  both  the  selection  rule  and  Bird’s  DSMC  method  show 
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excellent  agreement.  This  serves  to  establish  the  validity  of  the  selection  rule,  as  the  DSMC 
method  is  known  to  compare  well  with  experimentally  determined  shock-wave  profiles.  In 
addition,  the  equation  on  which  the  selection  rule  is  based  is  shown  to  be  directly  related 
to  the  time-counter  procedure  in  the  DSMC  method,  further  establishing  their  equivalence. 
Finally,  the  results  of  several  example  simulations  are  presented,  where  large  numbers  of 
particles  are  used  in  the  simulations,  to  demonstrate  the  greatly  improved  computational 
efficiency  of  the  present  method. 

II.  STEPS  IN  A  PARTICLE  METHOD 

In  outlining  the  steps  taken  in  a  particle  method,  it  is  most  convenient  to  consider  a 
region  of  space  free  of  boundaries  or  other  special  consideration.  Also,  we  assume  *he  entire 
space  has  been  conceptually  divided  into  a  suitable  arrangement  of  cells,  similar  to  the 
approach  taken  in  a  finite  element  method.  The  overall  procedure  in  a  particle  simulation, 
which  applies  at  each  time  step  At,  can  be  represented  1  _  the  following  six  steps. 

(i)  Advance  all  particles  in  space  at  cc  .stant  velocity  for  the  time  At  and  ignore  possible 
collisions. 

(ii)  Use  the  new  spatial  positions  and  sort  the  particles  into  individual  cells. 

(iii)  For  a  given  cell  in  physical  space,  choose  a  pair  of  particles  at  random. 

(iv)  Use  a  selection  rule  to  determine  whether  the  pair  formed  is  selected  for  collision.  The 
selection  rule  also  determines  the  total  number  of  sample  pairs  formed  for  each  cell  and 
the  distribution  of  colliding  pairs  over  relative  velocity.  In  order  to  simulate  a  real  flow, 
these  values  must  be  consistent  with  those  predicted  by  kinetic  theory. 

(v)  If  accepted  for  collision,  conserve  momentum  and  energy  while  colliding  the  pair  to 
find  their  new  states.  The  orientation  of  the  relative  velocity  vector  after  collision  is 
chosen  randomly  in  correspondence  with  the  scattering  properties  of  the  intermolecular 
potential  considered. 

(vi)  Return  to  step  (iii)  to  process  the  proper  number  of  pairs  in  each  cell  and  then  to 
process  all  cells. 


3 


The  above  listing  is  presented  primarily  for  purposes  of  discussion,  since  some  of  the  steps 
are  often  combined  or  altered  in  particular  applications. 

The  sorting  operation  in  step  (ii)  and  the  pair  selection  process  in  step  (iv)  prove 
to  be  the  most  troublesome  in  writing  efficient  code  for  a  computer  having  a  vector  or 
massively  parallel  SIMD  architecture.  The  time-counter  procedure  in  Bird’s  DSMC  method 
corresponds  to  step  (iv).  To  describe  it,  we  assume  the  particles  have  been  sorted  into 
individual  cells,  i.e..  steps  (i)  through  (iii)  have  been  completed.  The  time-counter  procedure 
consists  of  the  following  three  operations. 

(a)  Retain  the  randomly  chosen  pair  if  R  <  &Tg/(o'Tg)ma.x- 

(b)  For  each  pair  retained,  compute  r  =  2/(NncrTg). 

(c)  Collide  pairs  while  71  +  T2  +  T3  +  . . .  <  t  for  their  cell. 

In  step  (a),  R  represents  a  uniformly  distributed  random  number  in  the  range  [0,  1],  g  is  the 
relative  speed  of  the  particle  pair,  <rT  is  the  total  collision  cross-section,  and  max  denotes 
the  largest  value  for  the  cell.  In  step  (b),  n  is  the  number  density  and  N  is  the  number  of 
particles  in  a  cell.  In  step  (c),  the  subscripts  represent  all  the  pairs  accepted  for  collision  in 
a  given  cell,  and  t  is  the  current  time  measured  from  the  start  of  the  simulation. 

As  is  immediately  obvious,  the  procedure  contains  two  conditional  branching  state¬ 
ments,  one  in  (a)  and  one  in  (c).  Step  (c)  proves  to  be  more  difficult  to  vectorize,  because 
the  point  at  which  the  calculation  terminates  is  not  known  a  priori  and  it  may  be  widely  dif¬ 
ferent  for  different  cells.  In  the  application  of  the  algorithm,  one  discovers  when  to  terminate 
the  calculation  only  during  the  calculation  itself.  This  is  the  concept  of  data  dependency 
and  we  will  return  to  the  problem  it  presents  in  our  discussion  below.  With  regard  to  the 
physics,  it  is  clear  that  the  branching  logic  is  instrumental  in  determining  the  total  sample 
size  for  each  cell.  It  selects  the  fraction  that  is  to  be  retained  for  collision  and  it  selects  the 
distribution  of  relative  speeds  over  the  chosen  pairs  (i.e.,  large  values  of  aTg  generate  small 
values  of  r,  and  consequently,  more  of  these  are  needed  to  produce  the  required  sum). 

The  well-known  success  of  the  DSMC  method  in  being  able  to  predict  a  variety  of  ex¬ 
perimental  conditions  has  been  widely  recognized  as  providing  a  strong  argument  in  support 
of  its  use.  Its  theoretical  justification  is  based  on  the  argument  that,  when  the  above  steps 
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are  carried  out  a  large  number  of  times,  for  a  large  number  of  particles,  they  lead  to  the 
single-particle  collision  frequency  given  by 

0  =  nW^g,  (1) 

where  the  overbar  represents  a  mean  quantity  in  the  sense  discussed  by  Bird.4  Following 
similar  logic,  our  approach  is  to  return  to  the  basic  relation  from  which  the  well-known 
bimolecular  collision  rate  for  an  equilibrium  gas  is  derived  and  develop  a  relation  for  the 
nonequilibrium  state  that  leads  to  a  selection  rule  governing  collisions  that  does  not  contain 
the  data  dependency  represented  by  step  (c). 


III.  BIMOLECULAR  COLLISION  RATE 

The  bimolecular  collision  rate  is  usually  derived  for  the  case  of  an  equilibrium  gas  and 
a  hard-sphere  molecule.  This  is  because  the  velocity  distribution  function  is  known  for  the 
equilibrium  «tate  and  the  total  collision  cross-section  for  a  hard-sphere  molecule  is  clearly 
defined  and  finite,  and  the  associated  scattering  is  isotropic.  As  a  consequence,  the  algebraic 
steps  and  required  integrations  can  be  carried  out  in  a  straightforward  way.  We  will  assume 
a  hard-sphere  molecule  in  our  initial  development,  but  the  gas  will  not  be  assumed  to  be  in 
equilibrium. 

The  numerical  value  of  the  (integrated)  bimolecular  collision  rate  gives  the  total  number 
of  pair  collisions  occurring  between  two  species  per  unit  of  time  in  a  unit  volume.  However, 
in  order  to  handle  rarefied  flows,  specific  knowledge  concerning  the  functional  dependence 
of  the  pair  collision  rate  on  the  relative  velocity  is  needed  for  the  case  of  a  nonequilibrium 
gas.  To  obtain  this  information  one  must  start  with  the  basic  expression  for  the  collision 
rate  between  two  species  A  and  B  given  by5 

zABdCdZ  =  J^-fA(C)fB(Z)(7rd2AB)gdCdZ,  (2) 

where  nA  is  the  number  density,  fA  is  the  velocity  distribution  function,  and  C  is  the 
molecular  velocity,  all  for  specie  A.  Also,  dAB  is  the  average  diameter  for  the  two  molecules 
and  g  is  their  relative  speed.  For  convenience,  we  use  the  notation  dC  =  dCxdCydCz. 
The  Kronecker  delta  is  used  to  account  for  the  fact  that  AB  collisions  are  the  same  as  BA 
collisions  when  A  and  B  are  the  same  specie. 
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Equation  (2)  can  be  written  in  terms  of  the  relative- velocity  vector  (Z  —  C)  by  use  of 
the  transformation 

9  =  Z-C 

G  =  ( mAC  +  mBZ)/(mA  +  mB), 

where  g  is  the  relative- velocity  vector  for  the  colliding  pair,  and  G  is  the  center-of-mass 
velocity  vector.  For  this  case,  the  Jacobian  of  the  transformation  is  unity,  and  therefore,  we 
have  dC  dZ  =  dGdg.  On  introducing  the  transformation  into  (2)  and  integrating  over  G, 
to  obtain  an  expression  in  terms  of  g  alone,  one  encounters  the  following  integral  operator 

F(g)  =  r  Ia(G  -  £-g)  fB(G  +  —g)  dG ,  (3) 

J — oo  mA  rnB 

where  m*  =  mAmB/(mA  +  mB)  is  the  reduced  mass.  Therefore,  Eq.  (2)  can  be  written 

ZAB  d9  =  \+Tab  F^dg^d2*B)9,  (4) 

where  a  distinction  between  zAB  and  zAB  is  made  to  account  for  the  different  dimensions  in 
the  two  equations.  Equation  (4),  as  it  stands,  has  the  proper  form  to  develop  our  desired 
collision-selection  rule,  but  it  is  useful  to  first  review  the  fact  that  in  many  instances  F(g)  has 
a  very  simple  structure,  which  proves  to  greatly  improve  the  utility  of  applying  a  statistical 
method  to  (4). 

The  function  F(g)  describes  how  AB  pairs  are  distributed  as  a  function  of  their  relative 
velocity  vector  g ,  and  clearly,  it  is  a  function  of  three  independent  variables.  However,  we 
will  show  that  it  is  always  a  rather  smooth  function,  and  that  it  can  often  be  approximated 
by  a  function  of  the  relative  speed  g  alone,  i.e.,  a  single  independent  variable.  This  is  in 
contrast  to  the  velocity  distribution  function  /(C),  which  has  the  potential  for  being  a  rather 
complex  three-dimensional  function.  For  example,  it  assumes  a  bimodal  distribution  in  the 
upstream  portion  of  a  shock  wave.  The  purpose  of  reviewing  the  simple  structure  of  F(g) 
is  that,  in  using  a  statistical  sampling  scheme  to  evaluate  Eq.  (4),  a  far  smaller  number  of 
samples  can  be  used  than  would  otherwise  be  necessary,  and  this  clearly  represents  a  very 
important  practical  consideration. 


In  establishing  the  fact  that  F(g)  is  a  smooth  function,  let’s  first  consider  the  case  of  a 
single  specie  gas,  for  which  the  integral  corresponding  to  (3)  is  given  by 

F (9)  =  r  f(0  -  g/2 )f(G  +  g/2)  dG.  (5) 

J — OO 
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The  first  important  observation  we  make  concerning  (5)  is  that  it  produces  an  even  function 
of  the  relative  velocity  vector  g,  i.e.,  we  have 

F(g)  =  F(-g).  (6) 

We  note  that  this  is  a  general  result  and  it  does  not  depend  on  any  assumptions  whatsoever; 
it  holds  for  any  conceivable  velocity  distribution  function,  /. 

The  second  observation  is  that  the  integral  can  be  carried  out  exactly  for  the  case  of 
an  equilibrium  gas,  and  it  leads  to  the  expression 

In  this  case  F°(g)  is  a  function  of  the  scalar  magnitude  g  alone.  That  is,  the  function  is  not 
only  symmetric  but  it  is  also  spherically  symmetric. 

The  third  and  final  observation  we  make  concerning  the  integral  operator  (5)  is  that 
it  has  the  unexpected  property  of  satisfying  the  central  limit  theorem.  This  can  be  seen  as 
follows.  For  the  case  of  a  single  independent  variable,  the  convolution  integral  is  defined  by 

/*/  =  f/(*-o/m,  (8) 

J— 00 

and  a  known  result  of  the  central  limit  theorem  is  represented  by 

/*/*/*/*•••  — ♦  gaussian  . 

The  corresponding  operator  suggested  by  (5)  is  different  from  (8)  in  the  position  of  the  minus 
sign,  the  appearance  of  the  variable  x/2,  and  in  the  presence  of  the  argument  (£  +  x/2); 
therefore,  we  define  a  new  operator  by 

/#/  =  f°°  fit  ~  X/2 m  +  ,/2R.  (9) 

We  first  note  that  the  integrand  in  (8)  could  be  replaced  by  /(x/2  —  £)/(x/2  +  £)  and  the 
conclusions  of  the  central  limit  theorem  would  still  hold.  This  can  be  seen  by  a  simple 
transformation  of  variables  that  leads  back  to  Eq.  (8).  On  using  the  operator  defined  by  (9), 
we  can  now  write 

/#/#/#/#  •  •  •  =  (/#/)#(/#/)#(/#/)#  •  •  • 

=  e#e#e#e#-.. 

=  e*e*e*e*---  — ♦  gaussian  , 
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where  e  is  defined  to  be  an  even  function;  and  where  the  second  step  makes  use  of  (6), 
while  the  third  step  makes  use  of  the  fact  that  (8)  and  (9)  are  equivalent  when  /  is  an  even 
function. 

At  this  point  we  have  demonstrated  that  the  operator  defined  by  (5)  satisfies  the  central 
limit  theorem  for  the  case  of  a  single  independent  variable.  However,  it  is  easy  to  see  that 
the  results  can  be  generalized  in  the  following  way.  If  one  forms  the  Fourier  transform  of 
both  sides  of  (9),  makes  use  of  the  transformation  u  =  ((  —  x/2)  and  v  =  (£  +  x/2)  and 
notes  the  separation  of  variables,  one  obtains  the  relation  F(u)  =  where  u> 

is  the  transform  variable.  For  the  case  of  three  independent  variables,  one  introduces  the 
three-dimensional  transform  variable  u>  and  performs  exactly  the  same  operation  on  (5), 
which  again  produces  a  product  of  Fourier  transforms  given  by  f(—u>)f(u>).  It  is  easy  to  see 
that  repeated  applications  of  the  operator  defined  by  (9),  or  (5),  lead  to  repeated  products 
of  the  same  Fourier  transforms.  This  is  the  property  one  uses  in  proving  the  central  limit 
theorem,  namely,  one  shows  that  the  resulting  expression,  suitably  normalized,  tends  toward 
the  Fourier  transform  of  the  gaussian  distribution.6  The  same  principle  can  be  used  to  show 
that  Eq.  (3)  also  has  this  property,  however,  its  convergence  to  a  gaussian  is  not  as  rapid 
because  of  the  potential  asymmetry  present,  i.e.,  one  obtains  /a{~ w)/b(w). 

We  now  see  that  when  the  integral  operator  in  (5),  or  (3),  is  applied  a  large  number  of 
times  it  produces  the  gaussian  or  Maxwellian  distribution.  The  significance  of  this  result  is 
based  on  the  fact  that  the  central  limit  theorem  is  a  very  powerful  theorem,  in  the  sense  that 
only  a  few  applications  are  required  to  obtain  a  very  close  approximation  to  the  gaussian 
distribution.  The  quick  convergence  to  a  gaussian  distribution  is  best  seen  by  considering 
several  simple  examples.  If  we  assume  a  one- dimensional  asymmetric  function  for  /,  as  shown 
in  Fig.  1(a),  then  a  single  application  of  the  operator  given  by  (9)  leads  to  the  symmetric 
function  F  seen  in  Fig.  1(b).  Here  we  see  that  the  symmetric  function  F  is  beginning  to 
look  like  the  gaussian  distribution. 

A  second  example,  which  makes  use  of  a  three-dimensional  representation  for  /,  is 
shown  in  Fig.  2.  The  integration  required  by  (5)  can  be  carried  out  analytically  in  three 
dimensions,  if  /  is  assumed  to  be  given  by  a  spherically  symmetric  Cauchy  distribution 
(l/(a  +  x2)).  Normalization  requires  that  the  Cauchy  distribution  be  truncated,  because 
the  required  integral  is  otherwise  divergent.  The  truncated  Cauchy  distribution  used,  and 
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the  corresponding  equilibrium  distribution  having  the  same  total  energy,  are  shown  in  Fig. 
2(a).  The  result  of  the  application  of  the  operator  (5)  on  both  functions,  shown  in  Fig.  2(b), 
makes  it  clear  that  the  action  of  the  central  limit  theorem  is  so  strong  that  the  two  results 
for  F(g)  are  very  nearly  the  same.  For  use  in  our  discussion  below,  we  also  show  in  Fig.  2(c) 
the  corresponding  distribution  functions  for  the  relative  speed  (magnitude),  which  we  define 
as  R(g). 

The  observations  made  above  can  now  be  used  to  introduce  an  excellent  approxima¬ 
tion  for  (5),  which  can  then  be  used  for  a  better  understanding  of  the  exact  relation  (4). 
Because  of  the  properties  of  the  central  limit  theorem,  we  know  that  F(g)  tends  towards 
the  gaussian  distribution;  and  according  to  (6)  it  is  guaranteed  to  be  symmetric.  We  can 
therefore  approximate  (5)  by  a  spherically  symmetric  function.  However,  this  spherically 
symmetric  function  is  not  the  equilibrium  function  (7)  when  the  gas  is  out  of  translational 
equilibrium  (see  Fig.  2(c)).  The  approximation  suggested  is  therefore  the  replacement  of  the 
symmetric  function  F(g ),  which  is  a  function  of  three  independent  variables,  by  a  spherically 
symmetric  function,  F(g ),  which  is  a  function  of  only  one  independent  variable.  Using  this 
approximation  and  noting  that  the  differential  element  dg  becomes  \ng2  dg  for  the  case  of 
spherical  symmetry,  the  expression  for  the  speed-dependent  bimolecular  collision  rate  for  a 
single  specie,  hard-sphere  gas  is  obtained  from  (4)  and  becomes 

ZAB  d9  =  Y  [4tt g2F(g )  (fy]  {ird2)  g. 

On  introducing  the  following  definition 

R(g)  dg  =  4 ng2F{g)  dg ,  (10) 

and  on  introducing  the  time  interval,  A t,  our  approximate  relation  becomes 

ZAB  dgAt  =  \^-R(g)dg^  \{vd2)gAt}.  (11) 

The  partitioning  shown  in  (11)  and  the  inclusion  of  the  time  step  At  are  both  introduced 
for  convenience  below.  The  exact  form  of  (11)  is  of  course  Eq.  (4),  which  we  repeat  using 
the  same  format  employed  in  (11). 

zab  dgAt  =  U^F{g)dg}  [{vd2AB)gAt]  (12) 
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To  emphasize  the  fact  that  (11)  is  exact  as  well  for  the  equilibrium  state,  we  can  use  it  to 
derive  the  known  bimolecular  collision  rate  for  a  single  specie  hard-sphere  gas  in  equilibrium. 
Defining  Zab  to  be  the  integrated  bimolecular  collision  rate,  and  dropping  the  At  from  (11), 
we  have 

n2 

Zab  =  y(*^2)  Jq  9R(a)dg. 

For  the  equilibrium  state,  we  use  (7)  and  (10)  to  evaluate  the  integral,  which  yields  the 
well-known  result 

ZAb  = -^=(*d2)C,  (13) 

where  the  mean  molecular  speed  for  an  equilibrium  gas,  C,  is  given  by 

C  =  \J%kTlnm  .  (14) 

Because  Eq.  (11)  is  conceptually  simpler  than  the  exact  expression  (12),  we  start  with 
its  discussion.  The  left-hand  side  of  Eq.  (11)  represents  the  number  of  binary  collisions  that 
occur  in  a  unit  volume  in  the  speed  range  g  to  g  +  dg ,  in  the  time  interval  At.  The  right-hand 
side  of  the  equation  can  be  viewed  several  different  ways,  depending  on  how  one  chooses  to 
group  it.  For  large  n,  the  quantity  n2/ 2  is  approximately  equal  to  the  number  of  pairs  that 
can  be  formed  out  of  n  particles  in  a  unit  volume.  The  function  R(g)  is  the  distribution  of 
pairs  over  the  relative  speed  g.  Therefore,  the  quantity  enclosed  by  braces  is  the  number  of 
pairs  found  in  a  unit  volume  in  the  speed  range  g  to  g  +  dg.  Now,  the  expression  in  square 
brackets  is  the  volume  swept  out  by  the  hard-sphere  particle  in  the  time  A t  as  it  moves  at 
speed  g.  If  this  small  volume  is  viewed  in  conjunction  with  the  larger  unit  volume  being 
considered,  then  its  numerical  value  can  be  interpreted  as  the  probability  that  the  particle 
pair  will  be  found  in  the  small  volume  (this  is  the  same  basic  concept  used  in  setting  up  Eq. 
(2)).  Using  this  interpretation,  we  can  write  Eq.  (11)  as  follows 

zab  dg&t  =  S  R{g)dg  Pa,  (15) 

where  S  is  the  sample  size  (for  a  unit  volume)  and  P$  is  the  probability  of  a  binary  collision. 
In  this  case  we  are  using  the  two  relations 

S  =  n2/ 2,  (16) 
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Ps  =  ( xd2)g&t . 


(17) 


The  statistical  sampling  scheme  defined  by  Eq.  (15)  is  to  sample  a  unit  volume  S  times, 
saving  only  the  pairs  falling  in  a  given  speed  range,  and  thci.  ave  only  the  pairs  for  which  a 
uniform  random  number  in  the  range  [0,  1]  is  less  than  the  fraction  Ps.  The  number  of  pairs 
retained  is  the  number  zab^Q  At. 

Following  the  grouping  used  in  (15),  we  can  write  the  exact  relation  (12)  in  the  same 
form  and  it  leads  to  the  expression  given  by 

ZAB  dgAt  =  S  F(g)dg  P3,  (18) 


where  in  this  case,  we  have 

q _  nAnB 
1  +6ab' 

and  where  Ps  is  again  given  by  (17).  The  purpose  of  deriving  Eq.  (15),  which  is  an  approx¬ 
imate  result,  along  with  the  exact  result  (18),  is  to  emphasize  the  fact  that  (18)  frequently 
has  the  behavior  of  (15).  This  is  important  because,  in  a  statistical  sampling  scheme,  far 
fewer  samples  are  needed  to  sample  a  smooth  one-dimensional  function  than  for  a  complex 
three-dimensional  function,  assuming  the  same  level  of  accuracy.  In  view  of  the  fact  that 
practical  considerations  often  limit  one  to  sample  sizes  for  a  single  cell  ranging  from  10  to 
100  sample  pairs,  for  each  time  step,  it  is  clear  that  such  a  small  sample  size  is  only  useful 
in  sampling  a  smooth  one-dimensional  function.  However,  knowing  that  (18)  often  has  the 
symmetry  of  (15),  one  can  operate  on  (18)  with  a  relatively  small  sample  size  and  still  be 
confident  of  obtaining  good  results.  However,  for  those  cases  where  F(g )  is  believed  to  be 
more  complex,  then  it  is  clear  that  a  larger  number  of  sample  pairs  is  needed  for  each  time 
step,  if  one  is  to  obtain  reliable  results.  This  situation  makes  clear  where  the  transition  from 
a  general  approach  characterizing  kinetic  theory  problems  to  a  more  restricted  approach 
typically  addressed  by  particle  methods  is  taken. 


IV.  SELECTION  RULE 

In  practice,  where  it  is  desirable  to  set  the  number  of  particles  in  a  unit  cell  to  as  large 
a  value  as  possible,  the  corresponding  number  of  possible  pairs  becomes  much  too  large.  For 
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example,  for  n  =  100  we  find  S  =  n2/ 2  =  5,000.  Because  it  is  not  practical  to  sample  a 
cell  5,000  times,  one  must  settle  for  a  considerably  smaller  sample  size.  If  we  reduce  the  size 
of  our  sample,  5,  then  it  is  clear  that  Ps  must  be  increased  by  a  corresponding  amount  in 
order  for  the  product  to  yield  the  same  rate.  The  approach  we  will  take  is  that  one  is  free 
to  specify  the  sample  size,  5,  provided  that  Ps  is  adjusted  to  leave  the  product  (15)  or  (18) 
unaltered.  Because  of  the  adjustment  that  can  be  made  to  the  value  of  Ps,  we  chose  to  view 
it  as  a  selection  probability  as  opposed  to  a  collision  probability.  The  concept  of  a  selection 
probability  applies  to  the  model  problem,  while  the  concept  of  a  collision  probability  applies 
to  the  original  physical  problem.  On  using  an  arbitrary  value  of  S ,  the  corresponding  value 
of  the  selection  probability  (17)  becomes 

p*  “  { ras} (20) 

In  the  work  that  follows,  we  will  assume  a  single  specie  gas  to  minimize  algebraic 
complexity.  However,  the  work  can  be  repeated  for  the  case  of  two  species  in  a  straightforward 
way.  Use  of  P3  in  a  computational  scheme  is  simplified  considerably,  if  one  introduces  the 
definition  of  the  mean-free  path  length  for  a  hard-sphere  molecule  and  replaces  the  physical 
constant  ird2  by  reference  values  of  number  density  and  mean-free  path  length.  This  can  be 
done  because  the  hard-sphere  relation 

A  =  -=J- — — , 

y/2  n(7rd2) 

clearly  holds  at  all  points  in  a  flow.  This  step  leads  to  a  more  useful  expression  given  by 


When  the  simulation  of  a  wind  tunnel  flow  is  carried  out,  freestream  conditions  may  be  used 
to  specify  the  required  reference  values. 

Because  the  values  of  S  and  Ps  can  be  adjusted,  it  is  appropriate  to  review  combinations 
that  lead  to  optimal  computational  efficiency.  Recognizing  that  certain  data  dependencies 
can  eliminate  the  possibility  of  exploiting  any  form  of  parallelism,7  we  are  interested  in 
selecting  an  algorithmic  form  that  is  free  of  these  dependencies.  An  unfortunate  example 
of  such  a  dependency  was  identified  above  in  reviewing  step  (c)  of  the  DSMC  time-counter 
procedure.  It  is  such  data  dependency  that  must  be  avoided  if  the  advantage  of  parallel 
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structure  is  to  be  gained.  Because  no  single  expression  adequately  covers  the  concept  needed 
here,  we  will  describe  the  algorithmic  form  having  the  desired  properties  as  being  open  loop, 
as  opposed  to  closed  loop  which  requires  an  input  from  the  state  of  the  system.  We  turn  our 
attention  first  to  a  specification  of  S. 

A.  Natural  sample  size 

The  principal  difficulty  in  writing  efficient  code  for  a  computer  having  a  vector  or  a 
massively  parallel  SIMD  architecture  is  introduced  by  the  data  dependency  associated  with 
step  (c).  This  difficulty  can  be  averted  when  using  selection  rule  (22)  by  employing  any 
sample  size  that  is  proportional  to  n,  which  becomes  a  natural  sample  size.  For  example, 
n/2  sample  pairs  out  of  a  total  of  n  particles  in  a  unit  cell  is  one  such  case.  Consequently, 
it  proves  to  be  useful  to  specify  the  sample  size  (for  a  unit  volume)  as  follows 

S  =  K(n/ 2),  (23) 

where  K  is  a  constant  to  b^  determined.  Using  this  value  for  the  sample  size,  we  then  have 
from  (22) 


The  advantage  introduced  by  (23)  can  best  be  understood  as  follows.  Assume  the  entire 
listing  of  particle  data  is  read  sequentially  while  cell  by  cell  pairings  are  made  of  the  particles 
as  they  occur  in  the  listing.  This  automatically  leads  to  n/2  pairs,  where  n  is  a  different 
number  for  each  cell.  For  K  =  2,  the  listing  is  processed  twice.  However,  in  order  to  assure 
random  parings  on  successive  reads,  a  random  entry  and  a  random  stride  (spacing  between 
entries)  is  actually  used.  Here,  the  number  of  pairs  created  is  tied  directly  to  the  population 
of  a  cell  and  is  a  number  that  automatically  arises  rather  than  a  number  imposed  on  the 
system.  Consequently,  selection  rule  (24)  may  be  applied  in  step  with  the  processing  of  the 
list  and  one  has  an  open-loop  algorithm.  We  also  note  that  in  the  application  of  (24)  it  is 
important  to  use  the  time  averaged  value  of  n/nre f,  to  reduce  statistical  fluctuations,  and 
therefore  it  is  a  known  quantity.  A  discussion  of  the  fact  that  K  need  not  be  an  integer  value 
and  how  that  case  is  handled,  together  with  details  on  a  Cray-2  implementation,  can  be  found 
in  a  thesis  by  McDonald.8  A  straightforward  application,  without  regard  to  optimization  of 
code  performance,  is  to  use  (24)  in  step  (iv)  to  select  the  pairs  that  collide. 
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Clearly  Eq.  (24)  can  only  be  used  as  a  probability  when  Ps  <  1,  and  therefore,  any 
application  must  account  for  this  restriction.  The  condition  Ps  <  1  can  now  be  used  to  fix 
the  minimum  permitted  value  of  K,  leading  to 


/  nma.x  \  5max  At 
V  nref  /  V2  Aref 


(25) 


In  actual  practice  the  value  of  Ps  is  monitored  during  a  simulation  to  see  that  it  remains  less 
than  unity  for  the  conditions  used.  However,  we  can  judge  these  conditions  by  estimating 
the  value  of  K,  which  requires  that  we  review  the  individual  quantities  appearing  in  (25).  A 
study  of  Eqs.  (7)  and  (10)  shows  that,  for  an  equilibrium  gets,  99.4%  of  the  pairs  fall  in  the 
range  0  <  g  <  gma.x,  where 

5m  ax  =  2.5^4  kT/m.  (26) 


On  frst  sight,  it  appears  that  (26)  is  unbounded,  and  consequently,  (25)  would  predict 
an  unacceptably  large  value  for  K  when  the  flow  Mach  number  is  large.  However,  the 
temperature  does  not  become  unbounded  in  a  simulation,  and  this  can  be  seen  as  follows.  An 
estimate  for  the  maximum  temperature  in  a  flow  can  be  found  by  considering  the  stagnation 
point  on  a  blunt,  adiabatic  body  in  an  Euler  flow.  Along  the  stagnation  streamline,  the 
energy  equation  leads  to  the  relation 

a* "  (nr + s? ) ul 


where  a3  is  the  speed  of  sound  at  the  stagnation  point,  M\  and  ui  are  the  upstream  Mach 
number  and  fluid  speed,  respectively,  and  7  is  the  ratio  of  specific  heats.  On  combining  the 
last  two  equations,  we  obtain  an  estimate  for  <7max, 

?ra"  =  5 +  Wf)  (27) 

This  last  equation  shows  that,  for  a  monatomic  gas  and  in  the  strong  shock-wave  limit,  we 
have  5max  =  Vo  ui.  Because  the  quantity  uj  is  normally  set  in  a  simulation,  it  is  easy  to  see 
that  <7max  roughly  twice  as  large  as  t»i,  and  therefore,  it  is  not  an  unbounded  quantity,  as 
one  would  initially  be  lead  to  believe  on  first  seeing  (26). 


Because  we  will  be  considering  shock-wave  profiles  in  a  monatomic  gas,  we  continue 
with  7  =  5/3  and,  in  the  strong  shock- wave  limit,  Eq.  (25)  yields  the  following  estimate 

K  «  6.3uiAt/Ai, 
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where  a  density  ratio  of  4  was  used  in  (25).  On  the  basis  of  our  experience  in  using  the  present 
method,  a  reasonable  value  for  the  dimensionless  group  corresponds  to  ujAt/Ai  =  0.15. 
This  leads  to  a  if  value  near  unity,  which  is  a  very  practical  value.  Because  the  hard-sphere 
model  leads  to  larger  K  values  than  power-law  models  for  realistic  gases,  we  will  defer  further 
discussion  of  this  topic  until  Section  V. 

B.  Constrained  selection  probability 

If  we  are  interested  in  a  more  automatic  procedure  for  limiting  the  selection  probability 
so  that  it  is  always  less  than  unity,  then  we  may  make  use  of  the  approach  introduced  by 
Bird,  that  is,  for  the  case  of  the  hard-sphere  molecule,  to  set 

—  9 1 5m ax  >  (28) 

where  gm*x  is  the  maximum  value  of  the  relative  speed  in  the  cell.  This  is  the  concept 
employed  in  step  (a)  of  Bird’s  time-counter  procedure.  With  this  choice,  one  is  assured  that 
the  selection  probability  never  exceeds  unity.  In  this  case,  the  sample  size  is  given  by 


On  comparing  the  expression  in  square  brackets  with  the  value  defined  by  (25),  it  is  easy 
to  see  that  the  two  are  directly  related,  as  they  must  in  order  to  satisfy  (18).  The  most 
important  fact  concerning  (29)  is  that  the  sample  size  is  imposed,  as  a  direct  consequence 
of  (28),  and  therefore  the  number  of  times  one  must  carry  out  step  (iv)  is  set  by  this  value, 
i.e.,  one  does  not  have  an  open-loop  algorithm.  Because  of  this,  the  creation  of  efficient  code 
for  a  computer  having  a  vector  or  a  massively  parallel  SIMD  architecture  becomes  rather 
difficult  when  this  algorithmic  approach  is  used.  Likewise,  in  the  application  of  (29),  one 
should  use  the  time  averaged  value  of  n/nre{  to  reduce  statistical  fluctuations. 

C.  Selection  rule  verification 

A  straightforward  test  of  the  theory  is  to  compare  computed  shock  wave  profiles  using 
the  present  approach  with  the  corresponding  profiles  predicted  by  Bird’s  DSMC  method.  In 
order  for  the  comparison  to  be  a  rigorous  test  of  the  selection  rule  alone,  free  of  uncontrolled 
effects,  we  replaced  the  time  counter  procedure  used  in  step  (iv),  when  using  the  DSMC 
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method,  by  selection  rule  (24)  and  followed  effectively  the  same  series  of  steps  (i)  through 
(vi)  in  carrying  out  the  simulations.  These  results  therefore  do  not  provide  comparative 
information  on  vector  versus  serial  code  run  times.  Fig.  3  presents  comparisons  of  shock- 
wave  profiles  for  Mach  numbers  of  3  and  10.  Both  the  density  and  temperature  profiles 
are  shown  in  each  case.  The  simulations  employing  the  present  method  were  carried  out 
using  selection  rule  (24),  with  K  «  1,  Ai  =  5  cell  lengths,  and  u i  A t  =  0.3  cell  lengths. 
The  upstream  average  particle  number  density  used  was  approximately  90  particles  per  cell 
and  the  results  were  time  averaged  over  1,000  time  steps.  The  DSMC  results  were  obtained 
from  shock-structure  data  kindly  provided  by  Boyd9  from  his  recent  work  at  NASA-Ames 
Research  Center.  It  is  clear  from  the  comparison  that  the  two  methods  give  nearly  the  same 
results.  This  has  been  found  to  hold,  as  well,  for  other  Mach  numbers  not  shown  here. 


V.  POWER-LAW  MOLECULES 

All  of  our  results  have  been  developed  for  the  case  of  a  hard-sphere  model  of  molecular 
interaction.  In  order  to  generalize  these  results,  we  must  return  to  Eq.  (2)  and  replace  the 
expression  for  the  total  collision  cross-section  for  the  hard  sphere,  given  by  n d\B,  by  the 
more  elementary  form  cr(g,  x)  dCl,  where  cr  is  the  differential  collision  cross-section,  \  is  the 
scattering  angle,  and  dfl  is  the  differential  solid  angle.  For  the  hard-sphere  model,  cr  is  a 
constant  and  the  scattering  is  isotropic.  For  the  more  general  case,  one  must  consider  the 
role  played  by  the  scattering  angle  as  well.  However,  in  the  application  of  a  particle  method, 
one  may  first  use  the  integrated  collision  cross-section,  crT(g),  in  (2)  and  then  choose  the 
orientation  of  the  relative  velocity  vector  after  collision  to  be  a  random  quantity,  selected  in 
correspondence  with  the  scattering  properties  of  the  molecular  interaction  considered.  This 
is  where  the  \  dependence  is  introduced,  which  also  explains  step  (v)  in  Section  II,  where 
the  procedural  steps  in  a  particle  method  were  outlined.  On  this  basis,  we  generalize  Eq.  (2) 
by  replacing  the  quantity  xd?AB  by  the  total  collision  cross-section  crT(g). 

Because  the  algebraic  steps  which  led  from  Eq.  (2)  to  Eq.  (12)  are  not  affected  by  the 
introduction  of  the  function  crT{g),  the  more  general  form  of  (12)  can  be  written  immediately 
as 


ZAB  dgAt  =  {  (rff~B)F^d9}  \-aA9)9&A- 


(30) 
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The  key  algebraic  step,  in  developing  the  dimensionless  form  (22),  was  the  introduction  of 
the  dimensionless  group 

G{g)  =  V2  nrefAre{  aT(g),  (31) 


which  was  obtained  from  (21)  for  the  case  of  a  hard-sphere  molecule.  In  anticipation  of  using 
the  same  algebraic  step  and  for  the  case  of  a  single  specie  gas,  we  write  Eq.  (30)  in  the  form 


zab  dgAt = { iF(a)  M(^) G(j) 


(32) 


The  algebraic  expression  for  Ps ,  for  the  open-loop  algorithm,  is  then  obtained  directly 
from  (32)  and  it  reads 

:)  (33) 


Pa  =  4  (—']  G(g)-J^~  . 
F  Vnref  /  y/2  Aref 


Likewise,  by  factoring  gcrT(g)  /  gm^xcrT(gmiLX)  from  (32),  the  sample  size  S  associated  with  the 
constrained  selection  probability  is  obtained  in  the  same  way  and  is  given  by 

.  5max  At 


S  =  H 
2 


nlf)  G{9max)  v^Aref 


(34) 


In  the  work  that  follows,  we  will  find  a  need  for  a  general  expression  for  the  mean-free 
path  length  in  a  single-specie  gas.  The  desired  relation  is  found  by  use  of  the  basic  definition 


0A  =  C, 


(35) 


where  0  is  the  single-particle  collision  frequency,  and  through  the  definition  of  the  collision 
frequency  expressed  in  terms  of  the  bimolecular  collision  rate,  given  by 


2  f°° 

0  =  ~  /  zab  dg. 
n  00 


The  factor  of  2  appears  because,  in  the  case  of  like  molecules,  each  collision  terminates  two 
free  paths.10  On  using  (30)  for  a  single  specie  gas,  the  single-particle  collision  frequency  is 
therefore  given  by 

/oo 

gvT(g)F(9)dg,  (36) 

•oo 


which  yields  the  required  expression  for  A  when  substituted  into  (35). 
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A  common  approach  in  generalizing  hard-sphere  results  is  to  assume  a  repulsive  force 
acting  between  molecules  that  is  produced  by  the  following  potential 

V  =  ar~a,  (37) 

where  r  is  the  distance  between  the  colliding  pair,  and  a  and  a  are  constants,  characterizing 
the  molecule.  The  principal  problem  that  arises  in  using  (37)  is  that  its  total  collision  cross- 
section  is  infinite  for  all  finite  values  of  a.  Because  of  this,  the  collision  frequency  and  the 
mean-free  path  length  are  not  clearly  defined  for  such  potentials. 

This  difficulty  has  been  addressed  in  several  different  ways:  1)  on  the  basis  of  certain 
physical  arguments,  a  cutoff  can  be  introduced  to  limit  the  range  of  the  potential  and  thus 
limit  the  size  of  the  computed  quantities;  2)  physical  arguments  can  be  used  to  replace  the 
total  collision  cross-section  by  either  the  total  cross-section  for  momentum  transfer,  aM,  or 
the  total  cross-section  that  arises  in  the  calculation  of  the  coefficient  of  viscosity,  crM;  or  3)  use 
of  a  novel  approach  introduced  by  Bird  in  which  the  isotropic  scattering  of  the  hard-sphere 
model  is  retained,  but  the  molecular  diameter  is  allowed  to  vary  as  a  function  of  g.  This 
model  is  called  the  variable  hard-sphere  (VHS)  model.  The  argument  put  forward  in  favor 
of  the  VHS  model  is  the  observation  that,  for  most  flows  of  interest,  the  variation  in  the 
collision  cross-section  has  a  far  greater  influence  on  the  structure  of  a  flow  than  any  variation 
in  the  molecular  scattering  characteristics. 

A.  Variable  hard-sphere  model 

The  variable  hard-sphere  model  was  introduced  by  Bird  in  1980,  as  a  practical  approach 
to  the  solution  of  engineering  problems.11  Isotropic  scattering  is  assumed,  like  the  hard-sphere 
model,  but  its  total  collision  cross-section  is  allowed  to  vary  with  g  as  follows 

VvHsid)  =  A/92w,  (38) 

where  A  and  u>  are  constants.  The  constant  u>  is  is  related  to  the  exponent  of  the  interaction 
law  (37)  by  the  relation 

u  =  2/a.  (39) 

In  the  development  of  his  model,  Bird  obtains  the  following  explicit  expression  for  the  total 
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collision  cross-section 


<rT{g)  =  <7re{[m*g2/2(2  -u)kTie{]  u. 


which  is  equation  (9)  in  his  paper,  and  the  following  expression  for  the  mean-free  path  length 

a  =  (r/rrfr/[(2-«rr(2-u)v?'n<7«f],  («) 

which  is  equation  (10)  is  his  paper.  These  two  equations  can  be  combined  to  obtain  the 
dimensionless  group  (31)  for  the  VHS  model,  which  is  given  by 


r(,  _  1  (2  kTte{y 

G(9)  T(2  —  a»)  \  m*g2  )  ’ 


where  we  note  that  m*  =  m/2  for  a  single  specie  gas. 

On  substituting  (42)  into  (33),  we  obtain  the  following  simple  relation  for  the  selection 
probability  for  the  VHS  model 


gAf 

V2  Aref  ’ 


where  the  quantity  D(u>)  is  given  by 


D(w)  =  (7r/2)u,/r(2  —  u;). 


For  the  case  of  the  hard-sphere  molecule,  a  =  oo,  or  u  =  0,  and  T(2)  =  1;  thus  D(0)  =  1  and 
therefore  Eq.  (43)  reproduces  the  hard-sphere  result  (24).  For  the  case  of  Maxwell  molecules, 
a  =  4,  or  u  =  1/2,  and  T(3/2)  =  y/x  /2;  consequently,  D(  1/2)  =  y/2  and  Eq.  (40)  becomes 


K  (n«() 


f^ref  At 

Aref 


This  last  relation  shows  that,  for  the  case  of  Maxwell  molecules,  the  bimolecular  collision 
rate  is  independent  of  the  relative  speed  g,  which  is  a  known  result.  We  also  note  that  on 
comparing  (24)  and  (45)  that  the  quantity  g/y/2  is  replaced  by  CK f,  which  is  a  much  smaller 
quantity  when  the  Mach  number  is  large.  Therefore,  the  condition  Ps  <  1  is  much  easier  to 
satisfy  for  Maxwell  molecules  than  for  hard-sphere  molecules. 
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B.  Cross-section  for  momentum  transfer 

An  alternative  to  using  the  constants  implied  by  (40)  is  to  assume  that  one  may  use 
the  collision  cross-section  for  momentum  transfer  given  by12 

/9  _  \  2/a 

<rM(g)  =  2?r  g~4/aA!(a),  (46) 

where  a  and  a  are  the  constants  appearing  in  (37)  and  the  quantity  Ai(a)  is  a  pure  number, 
for  which  tabulated  values  (where  v  =  a  +  1)  can  be  found  in  Chapman  and  Cowling.13 
The  collision  cross-section  for  momentum  transfer  reduces  to  the  hard-sphere  value  when 
a  =  oo,  and  therefore,  it  is  more  suitable  for  use  here  than  the  cross-section  associated  with 
the  coefficient  of  viscosity,  which  reduces  to  2/3  of  the  hard-sphere  value.  On  comparing  (46) 
and  (40),  it  is  clear  the  two  functional  forms  are  identical,  except  possibly  for  the  constants. 

In  order  to  compute  the  dimensionless  group  (31),  we  need  an  expression  for  the  mean- 
free  path  length  at  the  reference  conditions.  This  is  where  (35)  and  (36)  prove  useful.  On 
substituting  (46)  into  (36),  assuming  the  equilibrium  state  (7)  and  using  definition  (10),  we 
obtain  x  2 

nKtenl  =  ^(^faAl(a)J^(^y^r(2-2-).  (47) 

Equations  (46)  and  (47)  can  now  be  used  to  eliminate  the  collection  of  constants  grouped 

4 

within  the  braces  in  (47)  to  obtain  an  expression  relating  nref0ref  and  crM(g)ga .  Finally, 
on  using  (35)  to  introduce  the  mean-free  path  length,  and  (31)  to  obtain  the  dimensionless 
group  G(g),  we  obtain 


On  using  (39)  it  is  clear  that  (48)  and  (42)  are  identical,  and  therefore,  the  exact  same 
selection  rule  is  obtained  for  the  two  models,  i.e.,  Eq.  (43)  applies  to  both  cases. 


C.  Sample  size 


Because  the  two  models  predict  exactly  the  same  results  for  G(g)  and  Eq.  (43)  applies 
to  both  models,  the  sample  size  for  each  is  obtained  from  (43)  and  we  have 


( 0(u)  ( is^L. 

\  nref  /  \9mtLX  /  \/2  Aref 


(49) 
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In  order  to  study  this  equation,  we  need  an  estimate  for  the  quantity  CVef/^max •  This  can 
easily  be  obtained  from  the  definitions  for  C,  given  by  (14),  the  definition  for  the  Mach 
number,  and  by  use  of  Eq.  (27),  which  lead  to 


^ref 

<7max 


(50) 


A  better  estimate  for  the  ratio  nm!ix/nie{,  than  the  constant  used  in  Section  IV,  is  to  use  the 
stagnation-point  density,  as  a  function  of  the  upstream  Mach  number.  After  these  equations 
are  collected,  the  resulting  expression  for  K  becomes  sufficiently  lengthy  that  a  graphical 
display  of  the  results  becomes  more  appropriate  than  an  algebraic  display  of  the  resulting 
equation. 


Fig.  4  presents  the  functional  dependence  of  K  on  the  power-law  constant  u>,  for  fixed 
values  of  upstream  Mach  number  M\,  and  for  7  =  5/3.  In  evaluating  (49)  the  dimensionless 
quantity  At/Ai  was  assigned  the  value  0.3,  which  we  have  found  to  be  a  reasonable  value. 
Typical  values  for  u ;  for  physical  gases  fall  in  the  mid  range  of  the  upvalues  plotted.  For 
example,  u  ss  0.25  for  the  monatomic  gas  argon,  while  uj  ta  0.22  for  the  diatomic  gas 
nitrogen.  It  is  clear  from  the  figure  that  for  physical  gases  we  may  use  K  values  near  unity 
and  that  the  situation  even  improves  as  the  freestream  Mach  number  increases.  It  is  also 
clear  that  the  hard-sphere  model  (u  =  0)  represents  the  most  severe  case,  while  the  Maxwell 
molecule  (u>  =  0.5)  represents  the  least  severe  case. 

A  study  of  Eq.  (49),  and  the  plot  in  Fig.  4,  show  that  the  sample  size,  K,  is  more 
sensitive  to  the  density  than  to  the  temperature  in  the  flow.  Our  estimates  were  made  for 
the  stagnation  point  of  a  blunt  body  and  a  thermally  insulated  wall.  If  the  wall  were  a  cold 
wall,  then  the  temperature  of  the  gas  at  the  stagnation  point  would  be  low,  but  its  density 
would  be  high.  It  is  extremely  difficult  to  obtain  general  analytic  results  for  this  case,  but 
several  exploratory  runs  seem  to  indicate  that  the  net  effect  is  to  increase  the  value  of  K  by 
as  much  as  50%,  over  those  found  in  Fig.  4,  when  the  Mach  number  reaches  30.  However,  the 
figure  shows  that  K  =  2  would  also  handle  this  case  for  typical  physical  gases,  and  therefore, 
the  basic  approach  is  the  same. 


Equation  (43)  gives  the  selection  rule  for  the  open-loop  algorithm  for  either  the  VHS 
model  or  for  the  model  using  the  collision  cross-section  for  momentum  exchange.  We  found 
in  Section  IV  that  the  hard-sphere  limit  of  this  selection  rule  predicts  the  same  shock-wave 
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profiles  as  those  given  by  Bird’s  DSMC  method.  A  further  test  of  the  present  approach 
is  to  study  the  opposite  limit,  namely,  the  Maxwell  molecule  limit.  The  selection  rule  for 
the  open-loop  algorithm  has  already  been  derived  for  this  case,  and  it  is  given  by  Eq.  (45). 
Fig.  5  presents  comparisons  between  the  shock-wave  profiles  obtained  using  (45)  and  results 
from  the  DSMC  method,  for  the  case  of  Maxwell  molecules,  a  monatomic  gas,  and  for  the 
shock-wave  Mach  numbers  of  3  and  10.  In  our  simulations,  we  used  K  «  1,  Ai  =  1.5  cell 
lengths,  and  u\At  =  0.3  cell  lengths.  The  upstream  average  particle  number  density  used 
was  approximately  110  particles  per  cell  and  the  results  were  time  averaged  over  1,000  time 
steps.  It  is  quite  evident  from  the  comparisons  that  there  is  essentially  complete  agreement 
between  the  two  methods.  The  DSMC  results  were  again  obtained  by  making  use  of  shock- 
structure  data  provided  by  Boyd.9 

The  shock  wave  provides  a  test  for  a  nonequilibrium  gas,  while  an  appropriate  test  for 
an  equilibrium  gas  is  to  compare  the  collision  frequency,  or  mean  free  path  length,  in  a  single 
cell  with  the  predictions  of  kinetic  theory.  Equation  (41)  shows  that  A  varies  as  Tl!2n~l  for 
a  Maxwell  molecule  and  as  n-1  for  a  hard  sphere  molecule.  Because  the  former  represents 
the  more  rigorous  test,  we  show  in  Fig.  6  the  results  of  several  example  simulations  for  A 
versus  n  for  two  values  of  the  sound  speed,  differing  by  an  order  of  magnitude,  and  for  a 
range  of  densities  representative  of  many  simulations  of  interest.  It  can  be  seen  that  the 
comparison  is  quite  good,  which  further  supports  the  view  that  the  proposed  selection  rule 
is  valid. 

D.  Time-counter  method 

Because  the  shock-wave  profiles  predicted  by  the  present  theory  are  essentially  identical 
to  the  profiles  predicted  by  the  DSMC  method,  it  is  clear  that  the  two  approaches  must 
be  closely  related.  The  comparison  that  suggests  itself  next  is  to  consider  whether  Eq. 
(30),  which  is  the  central  equation  in  the  application  of  our  method,  can  be  related  to  the 
time-counter  procedure  in  the  DSMC  method. 

If  we  cancel  the  At  in  (30),  assume  a  single  specie  gas  and  integrate  over  g,  we  obtain 
the  relation 


Zab  =  Y  <  9*t(s)  >» 


(51) 
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where  Zab  is  the  integrated  bimolecular  collision  rate  and  the  operator  representing  the 
mean  quantity  is  defined  by 


,°o 

<9°t{9)>=  /  9°T(9)F{g)dg. 

J—oc 


If  V  is  defined  to  be  the  volume  of  a  small  element  in  physical  space,  then  the  quantity 

r  =  l/(VZAB),  (53) 

represents  the  average  time  between  binary  collisions  in  the  volume  V.  On  substituting  (51) 
into  (53)  and  defining  the  number  of  particles  in  the  small  volume  to  be  N  =  Vn,  we  have 

Nn  <  gcrT(g)  >  ^ 

The  average  number  of  binary  collisions,  Ncou,  occurring  in  the  volume  V  in  the  time  At  is 
given  by  Ncou  =  At/r  and  therefore  (54)  can  be  arranged  to  read 

*=i,~U<Xw>)-  (55) 

Turning  our  attention  to  the  time-counter  method  and  steps  (a),  (b)  and  (c),  the  sum 
represented  by  step  (c)  can  be  written  for  a  single  time  step,  At,  as  follows 


=  £( 
«=i  x 


Nng<rT{g ), 


where  Ncou  is  the  same  quantity  introduced  in  (55).  Now,  the  operator  representing  the 
mean  quantity  appearing  in  the  denominator  of  (55)  is  based  on  the  distribution  F(g),  as 
defined  by  (52).  However,  the  statistical  sampling  procedure  used  in  steps  (a),  (b)  and  (c) 
does  not  lead  to  a  mean  based  on  F(g),  because  the  filtering  action  of  step  (a)  modifies  the 
distribution  to  one  that  is  proportional  to  g<7T(g)F(g).  Identifying  this  modified  distribution 
by  the  symbol  Fc(g),  its  representation  in  normalized  form  is  therefore  given  by 

p  _  9<rT(g)F(g)  .... 


Fc(g)  = 


<  9<rT{g)  > 


On  this  basis,  the  sum  in  (56)  is  described  by  the  following  mean  quantity 


At  —  Ncou 


{NngaT(g) 
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where  the  mean  is  based  on  the  distribution  Fc(g)  and  the  operator  used  in  (58)  is  defined 
by 


/— L- 

\NngtTT(g)  /  F 


■r  £(-N^M)F'(9)i9- 


(59) 


Noting  the  cancellation  that  occurs  when  (57)  is  substituted  into  (59)  and  the  fact  that  F(g) 
is  a  normalized  distribution,  we  see  that  Eq.  (58)  reduces  to 


At  (jVn  <  gaT{g)  >)  ’ 


(60) 


which  is  exactly  the  same  as  Eq.  (55).  Therefore,  the  two  approaches  are  entirely  equivalent, 
which  explains  why  the  comparisons  of  shock-wave  profiles  resulted  in  such  close  matches. 

The  key  step  connecting  the  two  methods  is  the  relation  between  the  operators  appear¬ 
ing  in  (55)  and  (58),  i.e., 

_ i _ 

<9^M>F  \9<7t(9 ) 

This  is  a  general  result,  which  can  be  derived  on  the  basis  of  probability  theory  alone,  and 
it  does  not  depend  on  the  particular  function  gcrT(g)  used  here.  Because  the  mean  and 
reciprocal  operators  do  not  commute,  Eq.  (61)  is  needed  to  make  the  proper  transformation. 


VI.  DISCUSSION  AND  CONCLUSIONS 

The  key  equation  in  the  application  of  the  present  method  is  Eq.  (30).  It  gives  the 
number  of  AB  collisions  occurring  in  a  unit  volume,  in  a  specific  relative-velocity  range,  in 
the  time  At.  The  selection  rule  for  the  open-loop  algorithm  which  follows  from  this  equation, 
for  the  case  of  a  single  specie  gas  and  the  VHS  model,  or  the  model  using  the  collision 
cross-section  for  momentum  exchange,  is  given  by  Eq.  (43).  Additionally,  the  corresponding 
equation  for  the  appropriate  sample  size,  used  with  the  selection  rule,  is  given  by  Eq.  (49). 
The  basic  difference  between  the  present  method  and  the  DSMC  method  is  in  the  application 
of  step  (iv).  Briefly,  in  the  present  method,  the  list  of  particle  data  is  processed  K  times, 
automatically  creating  S  =  K(n/2)  random  sample  pairs,  where  n  is  a  different  number 
for  each  cell.  On  using  selection  probability  (43)  in  step  with  the  processing  of  the  list,  to 
determine  which  pair3  collide  in  each  cell,  the  resulting  collisions  then  satisfy  (30)  in  terms  of 
the  number  of  collisions  per  unit  volume  in  the  time  At  and  in  the  distribution  over  relative 
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velocity  g.  The  sole  interest  in  the  present  method  is  that  it  leads  to  an  open-loop  algorithm 
and  greater  computational  efficiency  on  certain  machines,  as  opposed  to  the  closed-loop 
sequence  of  steps  required  by  the  time  counter  method  described  by  the  three  operations 
(a)  through  (c).  The  discussion  following  Eq.  (24)  should  be  consulted  for  a  more  precise 
description  of  the  method. 

The  aim  of  our  discussion  has  been  to  first  establish  the  validity  of  the  proposed  selection 
rule,  and  not  on  a  comparison  of  execution  speeds  for  different  computer  codes.  The  extensive 
programming  details  used  in  implementing  the  present  theory  on  the  Cray-2  supercomputer 
is  presented  in  Ref.  8  and  will  not  be  reviewed  here.  Nevertheless,  the  presentation  of  several 
example  cases  representing  the  capability  of  the  method  are  appropriate,  in  order  to  establish 
that  the  method  does  in  fact  offer  a  computational  capability  which  has  not  existed  before. 

Figs.  7,  8  and  9  present  three  such  examples.  In  each  of  these  cases,  the  flow  was 
simulated  as  a  three-dimensional  flow,  even  when  the  intent  was  to  study  a  two-dimensional 
problem.  Also,  all  the  cases  represent  an  ideal  diatomic,  hard-sphere  gas  with  an  adiabatic, 
no-slip  boundary  condition  on  the  body  surface.  The  Knudsen  number  in  each  case  was 
chosen  to  be  reasonably  small  so  that  near-continuum  flow  would  be  realized.  The  sample 
size  K  ss  1  together  with  the  parameter  U\At  =  0.3  cell  lengths  were  used  in  all  three 
cases.  In  addition,  all  the  simulations  were  carried  out  using  cubical,  three-dimensional  unit 
cells  everywhere  in  space.  This  we  consider  to  be  a  very  important  component  in  producing 
code  that  executes  very  quickly.  All  the  simulations  were  carried  out  on  the  NAS  Cray-2  at 
NAS  A- Ames  Research  Center. 

Fig.  7  presents  the  results  of  a  simulation  for  a  Mach  6  flow  about  a  circular  cylinder. 
The  simulated  wind  tunnel  consisted  of  a  space  of  120  cells  high,  80  cells  long  and  3  cells 
deep,  and  the  circular  cylinder  was  40  cells  in  diameter.  The  Reynolds  number  based  on  the 
diameter  was  720.  The  simulation  was  carried  out  using  1.0  x  106  particles  and  Ai  =  0.5  cell 
lengths.  The  figure  shows  the  density  distribution  in  the  flow,  where  the  display  is  in  the  form 
of  an  infinite  fringe  interferogram.  The  interest  in  the  problem  was  in  determining  whether 
the  correct  shock-wave  detachment  distance  would  be  predicted  by  the  simulation.  This  was 
found  to  compare  very  favorably  with  experiment,  giving  confidence  in  the  programming 
and  the  simulation  method  used.  The  simulation  was  carried  out  using  a  full  cylinder  and 
therefore  the  symmetry  seen  in  the  flow  gives  further  confirmation  of  the  method  employed. 
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Fig.  8  represents  a  simulation  for  a  Mach  6  flow  about  a  10  degree  half-angle  wedge, 
where  the  velocity  vector  field  is  displayed.  The  dimensions  of  the  simulated  flow  were  90 
cells  high,  256  cells  long  and  3  cells  deep.  The  wedge  was  128  cells  long  and  22  cells  high; 
and  based  on  the  full  base  height,  the  Reynolds  number  was  1760.  In  this  case  a  total  of 
0.74  x  106  particles  were  used  in  the  simulation  with  Ai  =  0.23  cell  lengths.  The  aim  of  the 
investigation  was  to  see  how  detailed  the  resulting  velocity  field  would  be.  Of  great  interest 
is  the  detail  found  in  the  structure  of  the  boundary  layer  on  the  wedge  surface  and  in  the 
detail  found  in  the  region  immediately  behind  the  wedge.  The  vortex  flow  behind  the  wedge 
can  be  seen  very  clearly,  if  the  scale  of  the  flow  is  increased  to  allow  one  to  focus  on  it. 
This  result  contradicts  the  conclusion  drawn  by  Meiburg14  who  argued  that  the  molecular- 
dynamics  method  of  particle  simulation  must  be  employed  in  order  to  see  such  phenomena. 
Additional  flow  detail  which  can  be  seen  more  clearly  in  the  pressure  field  than  in  the  velocity 
field  is  the  existence  of  a  wake  shock. 

Fig.  9  presents  the  results  of  a  study  of  a  true  three-dimensional  simulation,  and  rep¬ 
resents  an  attempt  to  learn  how  well  a  three-dimensional  flow  can  be  resolved  with  the 
capability  at  hand.  The  simulated  wind  tunnel  consisted  of  a  space  of  120  cells  high,  60  cells 
long  and  60  cells  deep.  A  total  of  9.5  x  106  particles  were  used  in  the  simulation.  The  body 
diameter  was  44  cell  lengths,  Aj  =  1.0  cell  lengths,  and  M\  =  35.  The  Reynolds  number 
based  on  the  body  diameter  was  2300.  The  view  shown  is  the  pressure  field  in  the  central 
plane  of  a  blunt  nosed  body  inclined  at  an  angle  with  respect  to  the  freestream.  The  body 
represents  a  typical  blunt  shape  being  considered  for  use  on  spacecraft  for  aerodynamically 
assisted  orbital  changes.  These  results  were  obtained  as  part  of  a  study  at  NASA-Ames 
Research  Center  to  apply  the  methods  discussed  above  to  practical  aerospace  problems.15 

All  of  the  above  cases  were  run  at  a  program  execution  rate  of  around  1  to  2  xl0~6 
seconds  per  particle  per  time  step.  For  a  simulation  employing  106  particles,  approximately 
1  to  2  seconds  of  computation  time  is  required  per  time  step.  Because  a  statistical  method 
requires  time  averaging  over  many  time  steps,  to  reduce  statistical  fluctuations  in  the  results, 
the  total  run  time  depends  directly  on  the  length  of  time  averaging  employed;  thus,  the  use  of 
roughly  1,000  time  steps  leads  to  approximately  one-half  hour  of  Cray-2  computation  time. 
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of  Eq.  (5).  (a)  An  assumed  spherically  symmetric,  truncated 
anction  F(g)  .  (c)  The  function  R(g).  _ ,  test  function; 


(b) 


Fig.  3.  Shock-wave  density  and  temperature  profiles  for  a  monatomic,  hard-sphere  gas.  The 
symbols  designate  the  solution  obtained  using  selection  rule  (24)  and  the  solid  (density)  and 
dashed  (temperature)  curves  are  for  the  DSMC  method  (Ref.  11).  (a)  M\  =  3.  (b)  M\  =  10. 


Fig.  4.  Functional  dependence  of  the  sample  size  K  on  the  power-law  constant  u>  and  the 
freestream  Mach  number  M\,  for  the  case  7  =  5/3  and  uiAf/Ai  =  0.3.  M\  =  2  (solid  line), 
4,  8,  16,  32,  64  (dashed  line). 
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Fig.  5.  Shock- wave  density  and  temperature  profiles  for  a  monatomic,  Maxwell-molecule  gas. 
Symbols  designate  the  solution  obtained  using  Eq.  (45)  and  the  solid  (density)  and  dashed 
(temperature)  curves  are  for  the  DSMC  method  (Ref.  11).  (a)  Mi  =  3.  (b)  Mi  =  10. 
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Fig.  6.  Comparisons  of  the  mean-free  path  length  obtained  in  simulations  based  on  selection 
rule  (45)  for  Maxwell  molecules  and  with  kinetic  theory.  Sound  speed  (cell  lengths/time 
step):  o,  0.2;  +,  0.02. 


Fig.  7.  Flow  past  an  adiabatic  circular  cylinder  for  the  case  of  a  hard-sphere,  diatomic  gas 
and  a  Mach  number  of  6.  Density  distribution  is  shown  in  the  form  of  an  infinite  fringe 
interferogram.  A  total  of  1.0  x  106  particles  in  a  space  of  120  x  80  x  3  cells  were  used  in  the 
simulation. 


Fig.  8.  Velocity  vector  field  about  an  adiabatic,  10  degree  half-angle  wedge  for  the  case  of 
a  hard-sphere,  diatomic  gas  and  a  Mach  number  of  6.  A  total  of  0.74  x  106  particles  in  a 
space  of  90  x  256  x  3  cells  were  used  in  the  simulation. 


Fig.  9.  Pressure  distribution  in  the  central  plane  of  a  three-dimensional  Mach  35  flow  de¬ 
scribed  in  the  text.  A  total  of  9.5  x  106  particles  in  a  space  of  120  x  60  x  60  cells  were  used 
in  the  simulation.  The  body  diameter  was  44  cell  lengths  and  Ai  =  1.0  cell  lengths. 
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ABSTRACT 

The  flowfield  about  vehicles  re-entering  the  atmo¬ 
sphere  is  characterized  by  hypersonic  rarefied  flow 
in  thermochemical  non-equilibrium.  An  appropriate 
method  for  computing  such  flows  is  direct  particle  sim¬ 
ulation,  an  efficient  implementation  of  which  utilizes 
the  vector-processing  capability  of  current  supercom¬ 
puters.  The  fundamentals  of  extending  the  vectorized 
simulation  to  model  chemically  reactive  flows  are  de¬ 
scribed  here.  A  review  is  presented  of  the  reaction 
physics  of  dissociation,  exchange  reactions,  and  re¬ 
combination,  with  attention  to  internal  energy  mode 
contribution.  Reaction  selection  rules,  as  functions  of 
reactive  collision  energy,  are  reviewed.  Reaction  me¬ 
chanics  modeling  is  constrained  by  considerations  of 
detailed  balance  at  equilibrium,  and  conservation  of 
linear  momentum  and  energy.  The  proposed  means 
of  partitioning  post-reaction  energy  among  the  energy 
modes  of  the  products  is  based  on  proportional  energy 
exchange  per  mode.  Details  of  reaction  mechanics  per 
reaction  are  presented  with  particular  attention  to  the 
quantum  nature  of  the  vibrational  mode.  Verification 
of  the  models  is  made  through  simulation  of  super¬ 
heated  diatomic  reservoirs  relaxing  thermochemically 
to  equilibrium. 


NOMENCLATURE 

A,  B,  C, ... 

monatomic  chemical  species 

AB,AC, ... 

diatomic  chemical  species 

n 

number  density 

forward  rate  coefficient 

K 

equilibrium  concentration  coefficient 

k 

Boltzmann’s  constant 

T 

temperature 

reduced  mass  of  colliding  pair 

a 

exponent  in  power  law  potential 

<7 

constant  in  collision  probability 

9 

relative  speed  of  collision 

a 

constant  in  Arrhenius  expressions 

b 

exponent  in  Arrhenius  expressions 

X 

position  of  a  particle 

u 

velocity  of  a  particle 

G 

collision  center-of-mass  velocity 
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Ei  activation  energy  for  reaction 

e,  £,  £  energies 

P  probability 

/  energy  distribution  function 

c  number  of  degrees  of  freedom 

9V  vibrational  characteristic  temperature 

(3  constant  in  reaction  probability 

ti>  exponent  in  Pd  function 

<j>  exponent  in  Pr  function 

Ad  degree  of  dissociation 

%  uniform  random  number 

q,  i,  j  vibrational  quantum  levels 

Q  maximum  quantum  level 

Superscripts: 

post-collision  value 
post-dissociation  value 
equilibrium  value 

pair  of  species  a  and  b 
collision  pair,  A  and  B 
orbital-pair  for  A  and  B 
collision  partner  of  any  type 
collision 
dissociation 
recombination 

denotes  forward  reaction  rate 
denotes  equilibrium  concentration 
internal  energy  modes 
relative  translational  mode 
rotational  mode 
vibrational  mode 
reactions  in  equations  3  and  4 
energy  mode 
reverse  reaction 
reference  value 

INTRODUCTION 

Potential  development  of  vehicles  such  as  the  Na¬ 
tional  Aerospace  Plane  (NASP)  and  the  Aeroassist 
Flight  Experiment  (AFE)  have  renewed  interest  in 
modeling  hypersonic  rarefied  flow.  Such  flows  may 
be  characterized  by  non-equilibrium  between  molec¬ 
ular  energy  modes  and  between  the  concentrations  of 
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different  species  in  the  rarefied  gas.  Temperatures  ex¬ 
ceeding  20, 000/C  are  possible  behind  the  leading  shock 
of  a  vehicle  entering  the  atmosphere  with  near-orbital 
velocity  Such  conditions  excite  the  vibrational  mode 
of  molecules  in  the  air  leading  to  dissociation.  Ex¬ 
change  reactions  and  recombination  of  free  atoms  will 
occur  near  cold  catalytic  surfaces  and  during  expan¬ 
sion  of  the  gas  around  body  shoulders.  The  finite 
rates  of  thermal  excitation  and  reaction  processes  sig¬ 
nificantly  alter  the  character  of  the  flowfield,  affecting 
the  temperatures,  pressures,  and  heat-transfer  experi¬ 
enced  by  the  body,  and  thus  dramatically  influencing 
the  aerodynamics  of  the  vehicle. 

The  degree  of  rarefaction  of  these  flows,  at  alti¬ 
tudes  above  80km,  is  characterized  by  Knudsen  num¬ 
bers  exceeding  0.01.  The  low-density  nature  of  this 
flow  is  such  that  the  familiar  Navier-Stokes  equations 
are  not  applicable  due  to  failure  of  the  linear  consti¬ 
tutive  relations.  In  contrast,  discrete  particle  methods 
are  based  on  real-gas  molecular  models  and  provide 
potentially  powerful  alternatives  for  simulating  these 
flows.  One  frequently-applied  particle  simulation  is 
the  Direct  Simulation  Monte  Carlo  (DSMC)  method 
of  Bird  which,  in  a  recent  form1,  is  capable  of  solv¬ 
ing  only  modest-sized  problems  due  to  large  computa¬ 
tional  requirements.  An  alternate  particle  simulation 
method  developed  by  Baganoff  and  McDonald2-5  is 
structured  specifically  for  the  vector-processing  archi¬ 
tectures  of  current  supercomputers,  resulting  in  signifi¬ 
cant  improvement  in  performance.  Current  implemen¬ 
tation  of  that  simulation  models  the  three  dimensional 
non-reactive  flow  of  general  gas  mixtures4. 

This  paper  introduces  the  fundamentals  of  extend¬ 
ing  this  vectorized  particle  simulation  to  treat  chemi¬ 
cally  reacting  flows.  Comprehensive  modeling  of  chem¬ 
ical  kinetic  processes  in  a  gas  involves  solution  of  the 
Schrodinger  wave  equation,  and  is  therefore  ill-suited 
for  large-scale  simulations.  However,  meaningful  simu¬ 
lations  can  be  obtained  by  developing  reaction  models, 
applied  at  the  particle  collision  level,  which  yield  the 
correct  macroscopic  chemical  behavior.  The  types  of 
particles  used  in  the  simulation  and  their  properties 
are  presented  along  with  a  review  of  chemical  reac¬ 
tions.  The  basics  of  vector  processing  are  introduced 
with  a  brief  description  of  how  the  simulation  is  struc¬ 
tured  to  be  compatible  with  them,  including  discus¬ 
sion  of  the  mechanics  of  thermal  collisions.  The  reac¬ 
tion  selection  rules,  based  on  reaction  probabilities  as 
functions  of  collision  energy,  are  then  reviewed.  Re¬ 
action  mechanics  are  presented,  with  details  of  how 
post-reaction  energy  is  partitioned  among  the  energy 
modes  of  the  products.  The  thermochemical  model 
is  then  applied  to  a  reservoir  of  superheated  diatomic 
molecules  relaxing  thermally  and  chemically  to  equi¬ 
librium,  demonstrating  its  ability  to  simulate  reactive 
flows. 


REACTION  PHYSICS 

Atmospheric  air,  under  the  hypersonic  conditions 
described,  can  be  approximately  modeled  by  the  fol¬ 
lowing  five-species  set:  O?,  ^2,  NO,  O,  and  N .  A 
more  complete  model  would  include  some  of  the  ions 
such  as  NO+-  Associated  with  each  particle  in  the  flow 
is  its  position,  x  =  (r,  y,  z);  its  velocity,  u  =  (u,  v,  w)\ 
its  species-type;  and  the  energy  of  its  internal  modes, 
er<>t  and  tvn  if  it  is  diatomic.  Each  species  is  identi¬ 
fied  by  its  mass,  relative  diameter,  and  characteristric 
temperatures  of  vibration  and  dissociation. 

Associated  with  the  species  set  are  15  dissocia¬ 
tion/recombination  reactions,  and  4  exchange  reac¬ 
tions  suggested  by  Bird6  and  listed  in  Table  1.  In¬ 
cluded  in  the  table  are  the  parameters  pertaining  to 
the  macroscopic  reaction  rate  coefficients,  kf(T),  and 
the  species  concentration  coefficients,  K(T)7 ,  obtained 
by  fitting  experimental  data  to  Arrhenius  form.  Al¬ 
though  application  of  this  experimental  correlation  for 
temperatures  exceeding  10,000K  is  questionable,  it  is 
used  here  for  lack  of  a  better  alternative.  The  reaction 
selection  rules  Me  derived  to  reproduce  the  macro¬ 
scopic  temperature-dependence  given  by  the  Arrhe¬ 
nius  forms. 

Very  little  is  understood  regarding  the  non¬ 
equilibrium  behavior  of  a  chemically  reacting  gas. 
However,  reaction  rates  and  concentrations  are  well 
documented  for  gases  very  near  equilibrium.  There¬ 
fore,  to  model  non-equilibrium  flow  we  must  devise  re¬ 
action  models  based  upon  our  limited  knowledge  of  the 
microscopic  physics,  specifying  the  unknown  parame¬ 
ters  of  those  models  by  matching  the  experimentally- 
observed  macroscopic  equilibrium  behavior. 

Dissociation  reactions  proceed  as  indicated  by  the 
general  equation 

AB  +  X  —  A  +  B  +  X  (1) 

where  AB  is  diatomic,  A  and  B  are  monatomic,  and 
X  is  a  partner  of  any  type.  Dissociation  is  possible 
only  for  those  collisions  with  energy  exceeding  the  dis¬ 
sociation  threshold,  Ej.  This  represents  the  reaction 
energy  which  is  removed  from  the  energy  modes  con¬ 
tributing  to  the  reaction  and  stored  in  the  electron 
configuration  of  the  products. 

Exchange  reactions  are  typically  of  the  form 

AB  +  C->AC  +  B  (2) 

where  C  is  monatomic.  Again,  collisions  must  have 
energy  exceeding  the  threshold,  Ei,  to  initiate  the  re¬ 
action.  Associated  with  reaction  (2)  is  its  reverse  reac¬ 
tion  which  itself  has  a  threshold  energy,  Eir„.  During 
reaction  the  energy,  A  Ed  —  Ei  -  Eir„,  is  removed 
from  the  energy  modes  of  the  reactants  and  stored  as 
p  r.tial  energy  in  the  products. 

Recombination,  the  reverse  of  dissociation,  is  a 
three-body  process.  Benson  and  Fueno8  suggest  that 
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recombination  should  be  modeled  as  a  succession  of 
binary  collisions.  In  the  first  step,  given  by 

A  +  B-*(AB),  (3) 

two  atoms  collide  and  possibly  form  a  mutually  orbit¬ 
ing  pair,  denoted  in  parentheses  as  ( AB ).  The  pair  re¬ 
mains  intact  for  some  limited  lifetime  depending  upon 
the  collision  energy.  If  the  orbiting  pair  experiences 
a  collision  during  that  lifetime,  then  it  may  stabilize 
into  a  molecule,  completing  the  second  step  of  recom¬ 
bination,  given  by 

(AB)  +X  ^AB  +  X.  (4) 

Upon  stabilization,  the  reaction  energy,  Ed,  from  (1) 
is  removed  from  the  electron  configuration  of  the  re¬ 
actants  and  repartitioned  among  the  energy  modes  of 
the  products. 

If  no  stabilizing  collision  occurs,  the  orbital 
pair  splits  into  free  atoms,  representing  a  standard 
monatomic  thermal  collision.  Such  a  2-step  method 
has  been  adapted  for  this  simulation  as  described  be¬ 
low. 


SIMULATION  FUNDAMENTALS 

The  pipelining  architecture  of  a  vector  processing 
machine,  such  as  the  Cray  Y-MP,  permits  simultane¬ 
ous  operations  involving  collections  of  data,  yielding 
extremely  fast  computational  throughput.  However, 
algorithms  must  be  structured  accordingly  in  order  to 
utilize  this  capability.  Algorithms  which  involve  com¬ 
plex  conditional  branching  or  dependence  upon  other 
elements  within  a  single  vector  are  not  able  to  take 
advantage  of  this  architecture  effectively,  and  suffer  a 
loss  in  computational  throughput9.  When  processing 
a  list  of  data,  the  current  Stanford  simulation  handles 
conditional  branches  by  generating  new  subset  lists  for 
later  processing3.  The  positions,  velocities,  and  inter¬ 
nal  energies  of  each  particle  is  stored  in  the  particle 
array  in  computer  memory. 

The  basic  steps  in  a  particle  simulation  are  the 
collisionless  motion  of  particles,  followed  by  the  pair¬ 
ing  and  testing  of  neighboring  particles  for  possible 
collision.  The  collision  selection  rule  is  taken  from 
Baganoff  and  McDonald9-4,  and  is  based  on  the  power 
law  intermolecular  potential.  Each  class  of  interac¬ 
tion,  involving  species  a  and  6,  has  associated  with  it 
a  unique  power-law  exponent,  a„»,  and  cross-section 
constant,  <ra*,  as  used  in  the  collision  probability  of 
the  form 

Pc(9„6)  ~  <?a>9la;4/a"  (5) 


Particles  in  the  flow  are  paired  off  as  potential  colli¬ 
sion  candidates,  and  Pc  is  computed  per  pair  from  (5). 
To  select  colliding  pairs  from  among  all  candidates,  a 
uniformly  generated  random  number,  SR,  is  compared 
to  that  probability  such  that  the  collision  is  accepted 
if  3R  <  Pc.  Those  pairs  chosen  to  collide  are  placed 
in  a  new  list  which  is  later  processed  according  to  the 
collision  mechanics  of  the  method. 

An  algorithm  which  models  chemistry,  compatible 
with  this  list  generation  structure,  must  interrupt  the 
above  procedure  as  shown  in  figure  I.  After  pairs 
have  been  selected  for  collision,  they  must  be  tested 
for  reaction.  Depending  upon  the  class  of  collision, 
distinguished  by  the  types  of  species  involved,  collid¬ 
ing  pairs  are  divided  into  respective  class  lists.  Ideally 
each  class  should  then  be  tested  for  several  possible 
reaction  types  associated  with  it.  This  again  involves 
computation  of  probability  functions  as  discussed  be¬ 
low. 

Since  testing  all  possible  reactions  is  expensive  for 
a  given  colliding  pair,  a  simplification  is  employed  in 
this  simulation  which  significantly  reduces  the  amount 
of  conditional  branching  required.  That  is,  only  one 
type  of  reaction  is  considered  per  colliding  pair.  This 
simplification  is  appropriate  given  that  particles  are 
randomly  paired  and  the  sample  size  is  large.  Taking 
advantage  of  this,  the  type  of  reaction  considered  for 
the  pair  is  determined  not  only  by  the  types  of  particles 
involved,  but  also  by  the  specific  order  in  which  they 
appear  in  the  pair,  as  indicated  below: 


Dissociation  1 
Dissociation  2 
Exchange 
Form  Orbital  Pair 
Stabilize  Orbital  Pair 


AB  +  CD  —  A  +  B  +  CD 
AB  +  C  — *  A  +  B  +  C 
C  +  AB  — ►  AC  +  B 
A  +  B  -*  (AB) 

(AB)  +  X  —  AB  +  X 


After  reactive  collisions  have  been  selected,  they  are 
processed  according  to  the  reaction  mechanics  mod¬ 
els  appropriate  for  the  given  reaction-type.  Collisions 
which  do  not  react  are  processed  according  to  stan¬ 
dard  collision  mechanics  which  may  include  thermal 
relaxation  involving  possible  energy  exchange  among 
the  relative  translational,  rotational,  and  vibrational 
modes.  Such  collision  mechanics  are  discussed  in  de¬ 
tail  by  McDonald4  and  reviewed  briefly  below. 


THERMAL-COLLISION  MECHANICS 

A  collision  between  particles  a  and  b  has  associated 
with  it  a  center  of  mass  velocity,  G0*,  and  relative 
translational  velocity,  g„»  =  u„  -  u».  The  relative 
translational  energy,  t rei.t,  is  then  defined  by 


where  gat  is  the  relative  speed  of  collision.  The  power- 
law  exponent  typically  ranges  from  the  Maxwell- 
molecule  limit,  <*at  =  4,  to  the  hard-sphere  limit, 
where  <*„»  =  oo.  The  constant,  <t„»,  is  found  by  match¬ 
ing  some  specified  reference  mean-free-path4. 


«r «/.»  =  ^„i9at  (6) 

where  =  mam^/(ma  +  m»)  is  the  reduced  mass  of 
the  pair.  Neglecting  the  electronic  modes,  the  ener- 
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gies  available  in  a  given  collision  include  cre/.k  and  the 
rotational  and  vibrational  internal  energies  of  a  and  6. 

Mechanics  of  Elastic  Collisions 

After  an  elastic  collision,  the  particles  separate  with 
relative  translational  velocity,  g^4,  which  is  simply  a 
random  re-orientation  of  the  collision  approach  veloc- 
ity,  gai,.  This  represents  spherically  uniform  scattering 
of  post-collision  velocities  adapted  from  the  Variable 
Hard-Sphere  (VHS)  model  of  Bird10.  There  is  no  alter¬ 
ation  of  the  internal  modes  of  either  particle.  Note  the 
use  of  superscript  '  to  denote  post-collision  quantities. 

Mechanics  of  Rotational  Relaxation 

Some  collisions  may  permit  energy  exchange  be¬ 
tween  the  translational  mode  and  one  or  both  of  the 
internal  energy  modes,  promoting  relaxation  of  the  gas 
toward  thermal  equilibrium.  Those  collisions  involving 
the  rotational  mode  are  generally  more  frequent  than 
those  involving  vibration. 

Rotationally-relaxing  thermal  collisions  involve  re¬ 
partitioning  of  total  energy,  erei  +  (rot,  among  the  rel¬ 
ative  translational  and  rotational  modes  of  the  colli¬ 
sion.  This  is  adapted  from  the  method  of  Borgnakke 
and  Larsen11  by  sampling  the  fraction, 


directly  from  the  equilibrium  distribution ,  /(F). 
Knowing  the  types  of  particles  involved  in  the  collision, 
/(F)  is  easily  found  and  is  invariant  with  temperature 
when  assuming  that  the  translational  and  rotational 
energy  distributions  are  continuous.  The  post-collision 
rotational  energy  is  then  given  by 

('rot  =  ('r'l  +  (rot)  F.  (8) 

If  both  particles  in  the  collision  are  diatomic,  then 
('rot  is  randomly  divided  between  their  respective  ro¬ 
tational  modes.  A  probability,  Prot,  is  used  to  select 
which  collisions  will  include  rotational  relaxation. 

Mechanics  of  Vibrational  Relaxation 

A  probability,  Puitf  is  used  to  select  which  colli¬ 
sions  will  include  vibrational  relaxation.  These  colli¬ 
sions  are  modeled  by  McDonald4  in  a  manner  which 
avoids  direct  sampling  from  equilibrium  distributions 
such  as  /(F)  since  they  are  temperature-dependent 
when  vibrational  energy  is  involved.  That  model  will 
be  adapted  here  to  facilitate  its  use  in  reaction  me¬ 
chanics,  and  is  based  on  an  efficient  iteration  scheme. 
First,  energy  is  exchanged  between  rotation  and  vi¬ 
bration,  followed  by  an  exchange  between  rotation  and 
translation  via  rotational-relaxation  as  outlined  above. 
These  steps  may  then  be  repeated,  iterating  toward 
equilibrium. 

Note  that  consideration  must  be  made  of  the  quan¬ 
tum  nature  of  the  vibrational  mode.  For  the  simple 


harmonic  oscillator  we  define  vibrational  energy,  mea¬ 
sured  relative  to  the  ground  state,  ^ k6v ,  by 

ft iib  —  ?vi*(F0u)  (9) 

where  qvn  represents  the  number  of  quanta  associated 
with  a  given  vibrator,  and  k6v  represents  the  charac¬ 
teristic  vibrational  energy  per  quantum  level. 

To  model  the  random  mixing  of  rotational  and  vi¬ 
brational  energies  in  a  manner  which  promotes  equilib¬ 
rium,  we  take  advantage  of  the  fundamental  physical 
assumption  of  statistical  mechanics7.  This  states  that 
all  possible  divisions  of  internal  energy, 

(int  —  ( rot  A-  ( vii ,  (10) 

among  the  internal  modes,  ('rot  and  cFj,  are  a  priori 
equally  probable. 

First,  fjnt  is  quantized  by  the  characteristic  vibra¬ 
tional  energy,  k0v,  resulting  in  Q  +  1  quantum  levels: 
0,  1,  ...  Q.  The  limit,  Q,  is  found  from 


where  the  brackets,  J”,  denote  truncation. 

Upon  division,  the  outcome  quantum  level  for  vibra¬ 
tion,  q'vi j,  must  be  equally  likely  among  these  Q  -f  1 
quantum  levels.  A  means  of  doing  this  is  simply  to 
generate  a  floating-point  random  number  in  the  range 
[0,Q  +  1],  and  truncate  to  the  nearest  level,  as  given 
by 

[*(<5+ i)J  (i2) 

where  S  is  a  uniformly  generated  random  number  in 
the  range  [0, 1].  From  (12),  the  new  vibrational  energy 
given  by 

£v«*  =  Vvibi^v)-  (13) 

The  remaining  internal  energy  is  placed  in  the  rota¬ 
tional  mode  according  to 

(rot  £tnf  “  £wi6*  (14) 

As  proven  in  the  Appendix,  applying  this  method 
of  energy  partition  among  the  internal  modes  of  a 
diatomic  molecule  will  satisfy  detailed  balance  be¬ 
tween  the  vibrational  and  rotational  modes  at  equi¬ 
librium.  However,  when  out  of  equilibrium,  the  addi¬ 
tional  rotational-relaxation  steps  are  required  to  pro¬ 
mote  equilibrium  distributions. 

ENERGY  MODES  IN  REACTION 

For  dissociation  reaction  (1),  only  collisions  with  to¬ 
tal  energy  exceeding  the  activation  energy,  Ej,  are  ca¬ 
pable  of  dissociating  molecule  AB.  That  total  collision 
energy,  ee  =  Ecm,  is  the  sum  of  energy  from  all  con¬ 
tributing  modes,  m,  and  is  distributed  in  a  Boltzmann 
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distribution  with  Cc  fractional  degre  s  of  freedom12, 
given  by 


M*c)dec 


1 

r(Ce/2) 


(s/2  —  1 


e 


-<c/kT^££ 

kT 


(15) 


There  has  been  considerable  debate,  however, 
regarding  the  role  of  internal  energy  modes  in 
dissociation13.  It  has  been  unclear  which  of  the  en¬ 
ergy  modes,  available  in  a  given  collision,  contribute 
to  ec,  leading  to  confusion  over  the  appropriate  value 
of  Cc  Various  theories  on  chemical  kinetics  suggest 
that  reactions  occur  as  a  result  of  intimate  interac¬ 
tions  among  the  internal  energy  modes  of  the  reac¬ 
tants.  Sharma,  Huo,  and  Park14  suggest  that  disso¬ 
ciation  is  the  result  of  vibration-vibration  energy  ex¬ 
change  between  molecules.  Here,  AB  acquires  vibra¬ 
tional  energy  from  partner  X  until  it  reaches  its  reac¬ 
tion  threshold  and  dissociates.  Hansen15  suggests  *hat 
the  rotational  modes  contribute  to  vibrational  excita¬ 
tion  and  possible  dissociation.  Such  arguments  imply 
that  all  modes  contribute,  in  some  way,  to  fe.  In  addi¬ 
tion,  since  only  those  collisions  for  which  ee  >  Ed  are 
dissociation  candidates,  it  is  advantageous  to  include 
all  available  energy  modes  in  the  ec  sum,  resulting  in  a 
greater  number  of  sufficiently  energetic  collisions,  and 
thus  increasing  the  sample  size  for  reaction  selection. 

This  debate  concerning  internal  energy  mode  con¬ 
tribution  to  reaction  may  be  partly  resolved  in  a  self- 
consistent  manner  by  consideration  of  detailed  balance 
ar  ierived  by  Haas16.  More  involved  reaction  selection 
models,  fully  coupling  internal  modes  to  dissociation, 
have  been  studied  extensively1417  but  have  yet  to  ap¬ 
pear  in  large-scale  particle  simulations. 

Including  all  available  energy  modes  in  the  collision 
energy, 

€e  =  'rclAB.x  +  eintAB  +  Cintxi  (16) 
the  degrees  of  freedom  for  /e(fe)  is  simply  the  sum 
of  degrees  of  freedom  from  each  of  these  statistically 
independent  modes16, 

Ce=(4_— -  )  +  Ciw*4B  +  Cin«x>  (IT) 

V  aAB,x/rel 


where  the  subscripts  denote  the  energy  modes  con¬ 
tributing  to  C«-  In  (17)  the  4  -4/aAB  x  term  accounts 
for  the  relative  translational  energy  biased  by  the  col¬ 
lision  selection  rule  (5).  Due  to  the  rotational  and 
vibrational  internal  energies,  the  internal  degrees  of 
freedom  of  AB  is 


C»n(„a  —  ^  +  CvUab-  (16) 

Similarly,  C«nt*  18  zero  *f  X 's  monatomic,  and  is  given 
by  (18)  if  diatomic.  Normally,  the  vibrational  mode  at 
equilibrium  has  associated  with  it  Cvi*  degrees  of  free¬ 
dom,  defined  for  the  simple  harmonic  oscillator  by7 

_  2  ev/T 

('vti  ~  e#./T  _  1  •  l1®) 


Given  that  all  available  energy  modes  contribute  to 
the  reaction,  in  simulating  the  chemistry,  all  modes 
must  be  used  in  the  reaction  selection  rules,  and  all 
must  be  affected  by  the  reaction  mechanics,  as  pre¬ 
sented  below. 

DISSOCIATION  SELECTION 

Chemistry  modeling  involves  determining  which 
collisions  react  (reaction  selection;  and  then  perform¬ 
ing  the  reactions  (reaction  mechanics).  The  probabil¬ 
ity,  Pd,  that  a  given  collision  between  particles  AB  and 
X  results  in  dissociation  of  AB  as  in  (1)  is  computed 
as  a  function  of  the  energies  contributing  to  the  col¬ 
lision.  A  uniformly  generated  random  number,  SR,  is 
compared  to  that  probability  such  that  dissociation  is 
accepted  if  3R  <  Pd- 

The  expression  for  Pd  must  be  determined  such 
that  the  resultant  forward  reaction  rates  at  equilib¬ 
rium  match  those  from  experiment.  Possible  forms  of 
this  function  have  been  suggested15  based  upon  the 
physics  of  the  chemical  interaction  kinetics  and  are 
adapted  here  in  the  format  given  by  Bird18,  as 

*,(<.>  oo) 

By  matching  the  temperature  dependence  of  the 

simulated  dissociation  rate  due  to  application  of  (20), 
against  the  the  Arrhenius  fit19  for  the  rate  coefficient 
from  experiment,  given  by 

kJ(T)  =  afTtte-E*'kT  (21) 

the  free  parameters  in  Pd(fe)  are  found  to  be16 

(22> 

and 

.  (23) 

Vab.x  r(2~  r(^)  ** 

Note  that  ft  =  fiABiX  is  the  reduced  mass,  k  is  Boltz¬ 
mann’s  constant,  and  a,  and  bf  are  unique  constants 
per  reaction. 

With  the  exception  of  Cc,  all  parameters  in  Pd(ec) 
are  constant.  However,  and  thus  Cim  and  Cc  are 
very  weak  functions  of  temperature  for  the  range  of 
interest  in  hypersonic  flow,  as  shown  in  figure  2.  We 
therefore  fix  Cv<»  at  some  reference  temperature,  Tre/ , 
for  use  in  Pd-  An  appropriate  Tre /  is  any  intermediate 
translational  temperature  taken  from  that  portion  of 
the  flowfield  dominated  by  reactions.  Figure  3  shows 
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the  form  of  /«(««)  at  several  temperatures  along  with 
an  example  of  Pi(ec). 

Exchange  reactions  are  also  selected  via  Pi  with  the 
exception  that  Q„tx  *s  zero- 

RECOMBINATION  SELECTION 

As  noted  earlier,  recombination  is  modeled  by  a  2- 
step  process.  All  collisions  between  atoms  result  in  the 
formation  of  orbital  pairs  as  given  by  (3).  Even  for  a 
highly  dissociated  flowfield,  the  concentration  of  (AB) 
remains  very  low.  During  the  next  timestep,  all  {AB) 
which  collide  with  partner  X  are  tested  for  stabiliza¬ 
tion  by  the  recombination  probability  function 

Pr{(in,.,ein,x)  =  0r(tint.)-*°(cintxy*‘  (24) 

where  £;n(<>  and  (intx  are  the  internal  energies  of  the 
orbital  pair  (AB)  and  partner  X,  and  /?r,  <p„,  and  < px 
are  unique  constants  per  reaction.  For  monatomic  X, 
(intx  is  excluded  from  (24).  Note  the  use  of  subscript 
aon  to  denote  orbital  pairs,  (AB).  Those  collisions  se¬ 
lected  (by  comparison  of  a  random  number,  3i,  to  Pr), 
will  complete  the  recombination  process  of  equation 

(4)- 

The  form  of  Pr  in  (24)  was  chosen  because  it  leads  to 
the  correct  temperature  dependence  of  the  equilibrium 
concentration  coefficient,  K,  defined  for  reaction  (1) 
by 

K  =  (25) 

nAB 

where  K(T)  typically  assumes  an  Arrhenius  tempera- 
ture  fit  from  experiment, 

K(T)  =  aeTb‘e~B*lkT .  (26) 

The  parameters  in  Pr,  as  well  as  the  unknown 
interaction  potential  for  the  collisions  involving  or¬ 
bital  pairs,  a(AB)  x,  are  found  in  a  manner  similar  to 
that  employed  for  dissociation.  However,  a  complete 
derivation  involves  consideration  of  detailed  balance 
among  the  energy  modes  contributing  to  the  reaction, 


and  is  given  by  Haas16.  From  that  derivation, 
combination  parameters  are  found  to  be 

the  re- 

4 

°(AB).X  4-4/a 

*I<*AB.X  2V> - 

(27) 

4,-2  2 

(28) 

SC 

(29) 

and 

_  »(1 +*„..)  At(a,/at)4~^ 

F(^tjl) 

(30) 

r(2-£)r(2-£-*0)r(^-<M 

where  subscripts  3  and  4  correspond  to  the  reactions  as 
given  by  equations  (3)  and  (4),  respectively,  At  is  the 
duration  of  the  timestep,  hAB  is  the  Kronecker  delta, 
and  ip  is  given  by  (22). 

These  reaction  selection  rules  are  functions  of  col¬ 
lision  energies,  and  should  therefore  be  applicable 
throughout  the  flowfield.  The  free  parameters  of  these 
functions  were  determined  from  known  equilibrium  ex¬ 
perimental  data. 

REACTION  MECHANICS 

Once  reactive  collisions  have  been  identified  via  re¬ 
action  selection  rules  (20)  and  (24),  the  mechanics  of 
those  reactions  must  be  performed. All  reactions  of  in¬ 
terest  in  hypersonics  involve  some  alteration  of  the  en¬ 
ergy  modes  contributing  to  the  reaction.  The  objective 
here  is  to  determine  how  to  remove  or  add  energy  to 
each  of  the  energy  modes  involved  in  the  reaction  just 
prior  to  splitting  or  forming  particles  into  products. 
The  chemical-kinetic  details  are  not  well  understood 
regarding  energy  partition  among  the  energy  modes  of 
particles  during  reaction.  However,  the  modeling  of  re¬ 
action  mechanics  must  adhere  to  certain  fundamental 
constraints,  including  detailed  balance  at  equilibrium, 
and  conservation  of  momentum  and  energy. 

Proportional  Energy-Exchange  Model 

Recall  the  energies  contributing  to  dissociative  col¬ 
lision  (1)  as  identified  by  (16).  The  dissociation  energy 
must  be  removed  from  the  collision  energy, 

=  Ce  -Ei  (31) 

—  ^retxB.X  4*  £r otAB  ^vHab  4“  ^intx  Ei  (32) 
=  (relAB.x  +  ^rotAB  +  +  f<"‘x  •  (33) 

Note  the  use  of  superscript  "  to  denote  post-reaction 
quantities. 

Though  it  is  unclear  how  to  redistribute  the  ener¬ 
gies  of  (32)  among  those  of  (33),  a  convenient  means 
by  which  to  account  for  the  alteration  of  each  energy 
mode,  m,  is  to  proportionally  remove  the  reaction  en¬ 
ergy,  Ed,  resulting  in  post-dissociation  mode  energies, 
fm>  8>ven  by 

C  =  =  £m(l-^).  (34) 

€c  \  €c  ' 

This  method  is  simple  to  implement  and  can  be  ap¬ 
plied  to  exchange  reactions  and  recombination  as  well. 
Note,  however,  that  in  exothermic  reactions,  Ei  will 
be  proportionally  added  to  each  contributing  energy 
mode  rather  than  removed. 

At  equilibrium,  detailed  balance  dictates  that  there 
exists  no  net  energy  transfer  among  energy  modes  of 
all  species  in  the  gas7.  As  a  consequence,  the  reaction 
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models  must  promote,  or  relax  to,  equilibrium  for  a 
gas  in  thermochemical  non-equilibrium.  If  a  relaxing 
ensemble  of  particles  reaches  steady-state,  yet  exhibits 
non-equilibrium  behavior  such  as  a  difference  in  mode 
temperatures,  then  the  reaction  model  is  in  error. 

Just  as  in  the  method  of  Borgnakke  and  Larsen11, 
partitioning  post-reaction  energies  by  sampling  di¬ 
rectly  from  equilibrium  distributions  would  naturally 
accomplish  this.  However,  that  requires  that  these  dis¬ 
tributions  be  available  at  all  temperatures  for  all  pos¬ 
sible  reactions.  Alternatively,  the  proportional  energy- 
exchange  mode!  above  avoids  this  difficulty  since  it  re¬ 
quires  no  access  to  equilibrium  distributions,  but  fal¬ 
ters  in  that  it  will  not  readily  promote  equilibrium  for 
a  flow  initially  in  thermal  non-equilibrium.  It  is  there¬ 
fore  necessary  to  add  thermally-relaxing  steps,  as  in 
(7)-(14),  to  this  algorithm  for  reaction  mechanics. 

Mechanics  of  Dissociation 

Dissociation  reaction  (1)  is  divided  into  two  steps. 
First,  the  dissociation  energy  is  removed  from  the  re¬ 
actants,  creating  temporary  products,  as  described  by 

AB  +  X  — ►  AB"  +  X" .  (35) 

The  products,  AB"  and  X" ,  separate  with  relative 
translational  energy,  t"el AB  x,  modified  according  to 
the  proportional  energy-exchange  model  of  (34).  As 
mentioned  earlier,  decomposing  a  relative  energy  into 
post-collision  velocities  involves  random  re-orientation 
of  the  post-reaction  relative  velocity  vector,  g"B  x ,  as 
found  from  c"e,AS  and  equation  (6). 

If  X  is  diatomic,  the  post-dissociation  states  of  its 
internal  energies  are  not  found  individually  by  appli¬ 
cation  of  (34),  since  consideration  must  be  made  of 
the  quantum  nature  of  the  vibrational  mode.  Rather, 
the  internal  energy  as  a  whole  is  depleted  according 
to  (34), 

(36) 

and  is  then  redistributed  between  c"0(jf  and  f”nx  ac¬ 
cording  to  the  steps  employed  for  vibrational  relax¬ 
ation  in  equations  (11)-(14).  If  X  is  monatomic,  (jntx 
will  remain  zero. 

Now  that  AB"  and  X"  have  separated,  the  second 
step  in  the  dissociation  reaction  involves  splitting  AB" 
into  free  atoms,  given  by 

AB"  — >  A  +  B.  (37) 

In  order  to  conserve  linear  momentum  and  energy,  all 
of  the  post-dissociation  internal  energy  of  AB"  must 
go  into  the  relative  translational  energy  for  the  sepa¬ 
ration  of  A  and  B.  That  is, 

B  =  (intAB  =  ((rotAB  ~ (38) 


The  mechanics  for  exchange  reactions  (2)  are 
similar  to  that  for  dissociation,  replacing  Ed  by 
=  Ed  -  Edr.„- 

Mechanics  of  Orbital  Pair  Formation 

Orbital  pair  formation  is  the  first  step  in  recom¬ 
bination  modeling  as  given  by  (3),  where  two  atoms, 
A  and  B,  join  to  form  a  single  temporary  orbital  pair, 
(AB).  Conservation  of  linear  momentum  dictates  that 
the  velocity  of  (AB)  must  be  the  same  as  the  center  of 
mass  velocity  of  the  A  +  B  collision,  G a.b-  Conserva¬ 
tion  of  energy  requires  that  the  relative  translational 
energy,  ereiA  B,  must  be  stored  internally  in  (AB), 

=  (39) 

If  (AB)  remains  uncollided  on  the  next  time  step, 
it  is  split  into  free  atoms  A  and  B.  Doing  so  simply 
requires  that  (int(AB)  he  converted  back  into  treiA  B- 

Mechanics  of  Orbital  Pair  Stabilization 

Orbital  pair  stabilization,  given  by  (4),  represents 
the  second  step  which  completes  the  recombination  re¬ 
action.  The  reaction  energy,  Ed,  must  be  absorbed  by 
the  energy  modes  of  the  products,  such  that  the  post- 
reaction  collision  energy  is  given  by 

fe  —  f'cl(AB).x  f«"«(AB)  "F  (intx  "F  Ed 
—  f"  4-  r"  4-  f" 

=  fre(^a.x  +  Cin«4B  "F  ein»*  (40) 

According  to  the  proportional  energy-exchange  con¬ 
cept,  the  relative  translational  energy  of  separation 
of  stabilized  molecule  AB  and  collision  partner  X  is 
found  by  adapting  (34), 

1  +  — j.  (41) 

The  internal  energies  of  X,  if  applicable,  are  found 
in  the  same  way  as  in  dissociation  from  equation  (36), 
with  the  exception  that  reaction  energy  is  added  rather 
than  removed. 

Likewise,  the  internal  energy  of  AB,  given  by 

eintAB  =  eint{AB)  =  «m«MB)(l  +  ~)  (42) 

is  distributed  among  its  rotational  and  vibrational 
modes  according  to  the  method  employed  for  vibra¬ 
tional  relaxation  in  equations  (11)-(14).  However,  as 
stated  earlier  and  as  proven  in  the  appendix,  these 
steps  must  be  augmented  by  rotational  relaxation  to 
promote  equilibrium  distributions  for  flow  initially  out 
of  equilibrium. 
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RESERVOIR  SIMULATION 

A  single-species,  adiabatic  reservoir  of  molecules, 
superheated  in  translation  and  rotation  only,  will  re¬ 
lax  to  thermochemical  equilibrium  through  reactions 
and  energy-mode  exchanges.  The  final  temperature, 
T*.  and  composition  can  be  determined  analytically20. 
The  relaxation  history  of  this  reservoir  is  similar  to 
that  behind  a  normal  shock  wave. 

In  the  external  flows  typically  associated  with  hy¬ 
personic  re-entry  vehicles,  dissociation  and  exchange 
reactions  dominate  the  flowfield  chemistry17.  How¬ 
ever,  the  objective  here  is  to  test  the  simulation  in  an 
environment  in  which  all  types  of  chemical  reactions 
(dissociation,  exchange,  and  recombination)  are  pro¬ 
nounced,  such  as  would  occur  at  equilibrium  if  the 
temperature  and  diatomic  concentrations  remained 
high.  Consequently,  to  simulate  thermochemical  re¬ 
laxation  to  equilibrium,  a  model  species,  A 2,  is  used 
which  is  similar  to  02  except  that  ae  =  0.10  rather 
than  1200  as  given  in  table  1.  This  effectively  increases 
the  recombination  rate,  allowing  greater  concentration 
of  diatomic  particles  at  high  equilibrium  temperature 
where  reactions  are  significant.  Such  a  gas  represents  a 
“worst-case  scenario”  in  which  recombination  has  pro¬ 
nounced  effects  to  better  demonstrate  the  capability  of 
the  reaction  models.  The  reactions  pertaining  to  A2 
are 


perimental  fit,  where7 


and 

K  =  2(2nA3+nA)T^L-.  (46) 

These  plots  indicate  that  the  simulation  does  capture 
the  correct  equilibrium  behavior  of  a  relaxing  gas  as 
predicted  by  the  experimental  fits  of  the  reaction  rates 
and  species  concentrations  over  a  large  temperature 
range. 

The  simulation  was  then  run  for  the  5-species  mix¬ 
ture  to  demonstrate  its  applicability  to  hypersonic  air 
chemistry.  The  gas,  composed  of  79%  N2  and  21%  O2, 
is  superheated  in  translation  and  rotation  to  15,000K, 
and  allowed  to  relax  thermochemically  toward  equilib¬ 
rium.  The  plot  of  species  concentration  relaxation  is 
given  by  figure  8.  Note  the  rapid  increase  in  the  con¬ 
centration  of  atomic  oxygen  as  02  quickly  dissociates. 
Likewise,  dissociation  of  fV2  results  in  production  of  ni¬ 
trogen  atoms,  followed  eventually  by  formation  of  NO 
via  recombination  and  exchange  reactions. 


A2  +  A2  —  A  +  A  +  Aj  (43) 

Ai  +  A  —  A  +  A  +  A.  (44) 

As  plotted  in  figure  4,  A2  was  relaxed  from  16.871K 
to  steady  state.  The  mode  temperatures  correspond¬ 
ing  to  translation,  rotation,  and  vibration  are  plot¬ 
ted  in  time,  and  converge  to  the  equilibrium  tempera¬ 
ture,  analytically  determined  to  be  T*  =  7,000K.  The 
rates  of  internal  mode  relaxation  were  controlled  by 
invariant  inelastic  collision  probabilities,  Prot  =  0.20 
and  Pvn  =  0.02,  adapted  from  Bird20.  Energy- 
dependent  forms  of  these  probabilities,  suggested  by 
Boyd17,  would  better  model  thermal  relaxation  in  the 
transient  portion  of  the  flow,  but  are  of  little  benefit 
here  in  studying  equilibrium  behavior.  The  simulation 
started  with  40,000  diatomic  particles  and  ran  for  1500 
transient  steps  and  1000  steady-state  sampling  steps. 

The  relaxation  simulation  was  then  repeated  for 
temperatures  in  the  range  5,000K<  T*  <15,OOOK, 
with  reference  temperature  T,t /  =  10, 000K.  For  each 
case,  the  resulting  equilibrium  temperature  per  mode 
reached  by  the  simulation  is  compared  to  the  analytic 
value  in  figure  5.  The  close  agreement  between  mode 
temperatures  indicates  that  detailed  balance  is  main¬ 
tained,  as  is  necessary  for  equilibrium.  The  simulated 
forward  rate  coefficient  per  reaction  is  plotted  against 
T'  in  figure  6,  and  is  compared  to  the  experimental 
curve  of  k^T").  Likewise,  figure  7  compares  the  de¬ 
gree  of  dissociation,  A<j,  for  these  runs  against  the  ex¬ 


CONCLUDING  REMARKS 

The  chemistry  algorithms  presented  here  are  com¬ 
patible  with  the  existing  vectorized  simulation  of 
Baganoff  and  McDonald2-5.  The  dissociation  selec¬ 
tion  models  adequately  yield  the  correct  forward  reac¬ 
tion  rates  of  chemical  reactions  in  a  test  gas  as  found 
experimentally  near  equilibrium.  The  recombination 
selection  models  also  yield  correct  reaction  rate  be¬ 
havior  as  indicated  by  the  close  agreement  with  ex¬ 
perimental  fits  of  equilibrium  species  concentrations 
and  the  resultant  steady-state  temperature.  The  reac¬ 
tion  mechanics  models,  based  on  proportional  energy- 
exchange,  conveniently  account  for  the  effects  of  re¬ 
action  on  participating  particles  in  a  manner  consis¬ 
tent  with  the  constraints  of  detailed  balance  and  con¬ 
servation  of  linear  momentum  and  energy.  However, 
essential  to  these  routines  is  the  inclusion  of  thermal- 
relaxation  steps  which  promote  thermal  equilibrium. 

No  direct  attempt  was  made  here  to  model  the  tran¬ 
sient  of  the  relaxing  bath;  comparisons  to  known  be¬ 
havior  were  only  made  at  equilibrium.  However,  the 
chemistry  models  are  based  on  molecular  interaction 
physics,  independent  of  the  equilibrium  assumption 
which  is  used  only  to  evaluate  unknown  parameters. 
Keeping  in  mind  the  fundamentals  as  outlined  here, 
one  can  devise  more  sophisticated  chemistry  models, 
possibly  involving  more  species  and  greater  coupling 
of  internal  modes  to  reaction,  which  adequately  simu¬ 
lates  the  thermochemical  phenomena  of  interest. 
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APPENDIX 

Proof  of  Equilibrium  Convergence 

The  algorithm  described  by  equations  (11 )— ( 14)  for 
re-distributing  internal  energy  represents  random  mix¬ 
ing  of  continuous  and  discrete  distributions  pertaining 
to  the  rotational  and  vibrational  energy  modes,  re¬ 
spectively.  The  objective  here  is  to  prove  that  this 
algorithm  promotes  equilibrium. 

The  equilibrium  distribution  of  rotational  energy, 
for  which  (rot  =  2,  can  be  rewritten  from  (15)  in  the 
form 

/'(£)  =  e~‘de,  (47) 

where  e  =  erot/kT  is  normalized  energy,  and  super¬ 
script  *  denotes  equilibrium.  Quantizing  this  distribu¬ 
tion  with  respect  to  the  normalized  characteristic  en¬ 
ergy  of  vibration,  £  =  kOv/kT,  results  in  the  discrete 
equilibrium  distribution  over  vibrational  quantum  lev¬ 
els,  «, 


(«'+i)r 

P'(0=  J  Me)* 

it 

=  t/(l  -w)  (48) 

where  w  =  e~ 1  <  1  is  a  convenient  grouping  of 
terms.  For  a  given  rotational  energy,  e,  we  associate 
a  quantum  value,  j  =  [e/cj,  which  is  also  distributed 
as  in  (48),  such  that  the  total  quanta  from  (11)  is 
Q  =  i  +  j.  The  normalized  surplus  energy  remaining 
after  truncation  is  defined  by 


(  =  e  -  je,  (49) 


and  is  distributed  according  to  /*(f)  upon  renormal¬ 
izing  over  the  quantum  level,  j,  for  which  d(  =  de , 


ruw  = 

0+i)r 

/  f(e)de 

3 * 


-  e-f) 
e~* 


(1-w) 


4- 


(50) 


Note  that  /*(()  is  independent  of  rotational  quantum 
level,  j. 

Baganoff21  arrived  at  the  outcome  distributions, 
P{i')  and  P(j'),  in  the  following  manner.  Let  i'  and 
j'  denote  the  outcome  quantum  levels  due  to  uniform, 
random  division  of  total  quanta,  Q.  The  conditional 
distribution  for  given  Q,  is  therefore  constant  over 
the  Q  +  1  possible  outcomes, 


1 

Q  +  1 


(51) 


The  joint  distribution  of  i'  and  Q  is  then 

P(i',Q)  =  P(i'  \Q)  P;(Q)  (52) 

where  the  distribution  of  Q  is  given  by 

p;(Q)  =  P*(0)P*(Q)  +  P-(l)P'(Q  -  l)  +  ... 

Q 

=  £p-(i)P-(Q~i)  (53) 

i=0 

The  resulting  distribution  over  «'  is  found  by  summing 
(52)  over  all  Q, 

P(i')  =  £  X>-(»)P*(<3  -  i)P(i'  I  Q) 

Q=i'  i=0 

oo  Q  /.  \2  O 

(1  —  uiy 

= (i  -“)2  jt,  “Q 

Q=i‘ 

OO 

=  (i-u,)V'£X 

i=0 

=  u/(l-w)  (54) 

Comparison  of  (54)  to  (48)  proves  that  the  algorithm 
maintains  the  discrete  equilibrium  distributions  for  vi¬ 
bration  and  rotation. 

The  rotational  energy  after  mixing  is 

e'  =  A+C  (55) 


where  j'  =  Q-i',  such  that  de'  =  d£  over  the  quantum 
interval,  j*.  The  conditional  distribution  of  outcomes 
for  e‘  is  simply  the  joint  distribution  of  j'  and  (,  given 
Q,  found  by  multiplying  (51)  and  (50), 

f(e'  |  Q)de '  =  P(j'  |  Q)/;(0# 
e~(d£ 

*(Q+l)(l-w)-  (  } 

Following  the  same  development  as  in  (54),  we  solve 
for  the  distribution  over  e‘  by  summing  the  joint  dis¬ 
tribution  over  all  Q, 


f(e')de'=  £;  f;P’(i)P*(g-jW  I  Q)de' 
Q-i'  j=o 

_  ^  (1  —  w)  ujQ  e~(-  d£ 

Q=j‘j=0 


Q+ 1 


=  (1  —  w)  e-*  d(  ^ 
Q-i' 


=  (1  —  w)  e~t  up'  d£  T.  ujk 
t=o 

=  e-V'*»dt 

=  e~t‘de' 


P(i'  I  Q)  = 


(57) 
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Equations  (54)  and  (57)  prove  that  the  mixing  algo¬ 
rithm  maintains  both  the  discrete  and  continuous  dis¬ 
tributions  at  equilibrium. 

Note  that  (53)  represents  a  convolution  of  the  dis¬ 
crete  distributions,  P(i)  and  P(j).  By  virtue  of  the 
central  limit  theorm,  repeated  convolution  and  divi¬ 
sion  by  Q  +  1  will  result  in  relaxation  of  these  dis¬ 
tributions  toward  equilibrium  if  they  are  initially  out 
of  equilibrium.  However,  the  continuous  distribution, 
/(e')>  wiH  only  reach  equilibrium  if  both  P(j')  and 
/((s)  equilibrate.  Unfortunately,  the  distribution  over 
£  experiences  no  convolution  in  this  algorithm,  and 
therefore  will  not  relax  toward  equilibrium  unless  ad¬ 
ditional  rotational  relaxation  steps  are  included. 

To  test  the  algorithm,  a  reservoir  was  initialized 
with  no  vibrational  energy  and  with  rotational  energy 
uniformly  distributed,  and  was  then  thermally  relaxed 
via  the  mixing  algorithm  described  here.  Two  test- 
cases  were  run,  one  which  included  the  rotational  re¬ 
laxation  steps,  and  one  which  did  not.  The  initial  and 
final  distributions  for  rotation  are  shown  in  figure  9, 
along  with  the  final  distribution  for  vibration.  Note 
that  both  the  discrete  and  the  continuous  distribu¬ 
tions  are  Boltzmann  for  the  algorithm  which  included 
rotational  relaxation.  However,  when  rotational  relax¬ 
ation  was  excluded,  the  continuous  distribution  within 
each  quantum  interval  did  not  relax,  but  remained  uni¬ 
formly  distributed.  These  tests  prove  that,  for  a  gas 
out  of  equilibrium,  the  mixing  algorithm  will  relax  the 
discrete  distributions  to  equilibrium,  but  additional 
rotational-relaxation  is  needed  to  achieve  equilibrium 
for  the  continuous  distribution. 
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Table  1.  Reaction  equations  and  data6  7,19  for  the  5-species  air  model. 


REACTION 

at 

*/ 

a. 

be 

Ed/k  (K) 

.Vj  +  N2  **  N  +  N  +  JVj 

6.17  x  10-9 

-1.60 

18 

0.0 

113,000 

-V,  +  Oj  •*  N  +  N  +  02 

6.17  x  lO"9 

-1.60 

18 

0.0 

113,000 

iV2  +  NO  **  N  +  N  +  NO 

6.17  x  10“9 

-1.60 

18 

0.0 

113,000 

N2  +  N  ~  N  +  N  +  N 

1.85  x  lO"8 

-1.60 

18 

0.0 

113,000 

n2  +  o**n  +  n  +  o 

1.85  x  1Q-* 

-1.60 

18 

0.0 

113,000 

02  +  N2  s=*  O  +  0  +  iVj 

4.58  x  10"u 

-1.00 

1200 

-0.5 

59, 500 

02  +  02  0  +  0  +  Oj 

4.58  x  10-11 

-1.00 

1200 

-0.5 

59, 500 

02  +  NO  m  o  A-  0  +  NO 

4.58  x  10“u 

-1.00 

1200 

-0.5 

59, 500 

o2  +  n**o  +  o  +  n 

1.37  x  10-l° 

-1.00 

1200 

-0.5 

59, 500 

02  +  0  ^0+0+0 

1.37  x  lO"10 

-1.00 

1200 

-0.5 

59, 500 

NO  +  N2  5=*  N  +  0  +  N2 

3.83  x  lO"13 

-0.50 

4 

0.0 

75,550 

NO  +  02  *-  N  +  0  +  02 

3.83  x  10"13 

-0.50 

4 

0.0 

75,550 

iV0  +  JVO  *=*  N  +  0  +  WO 

3.83  x  10- 13 

-0.50 

4 

0.0 

75, 550 

NO  +  N*=*N  +  0  +  N 

7.66  x  10~13 

-0.50 

4 

0.0 

75, 550 

NO  +  O  **  N  +  0  +  0 

7.66  x  lO"13 

-0.50 

4 

0.0 

75, 550 

N2+0  ->  N0  +  N 

5.30  x  lO"17 

0.10 

37, 500 

NO  +  N  —  N2  +  O 

2.02  x  10"17 

0.10 

0 

02  +  N  —  NO  +  O 

5.20  x  10“” 

1.29 

3,600 

NO  +  O  -*  02  +  N 

3.60  x  10-« 

1.29 

19,700 

Note:  kf  from  equation  (13)  is  in  units  of  m3/(molecule-s),  K  from  equation  (18)  is  in  units 
of  mole/cm3,  and  T  is  in  Kelvin. 


Fig.  1.  Adaptation  of  a  vector-compatible  collision  algorithm  to  include  chemistry. 
F  iction-type  identification,  reaction-selection,  and  reaction-mechanics  are  added. 
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Equilibrium  Temperature,  T*  [K] 

Fig.  2.  Temperature  variation  of  degrees  of  freedom 
in  vibration  internal  modes  (Cm),  diatomic- 

monatomic  collsions  (Ceov),  and  diatomic-diatomic 
collisions  (Ceoo)  when  0V<»  =  2.270A. 


Normalized  Colisioo  Energy,  ^ 

Fig.  3.  Reaction  probability  function, 
and  collision  energy  distributions,  fe(ce),  where  ee  = 
(c/E, i,  for  dissociation  via  O*  +  Oj  collisions  at 
5,  000A,  10.000A-,  and  15.000A. 


Fig.  4.  Temperature  per  mode  per  species  in  A? 
thermochemical  relaxation  from  10, 871 A  to  7,OOOA. 


Fig.  5.  Simulated  vs.  analytic  equilibrium  temper¬ 
ature  per  energy  mode,  T*  in  A*  relaxation. 


Fig.  8.  Simulated  forward  rate,  kf,  vs.  equilib¬ 
rium  temperature  compared  with  experimental  fit  in 
A?  relaxation. 


Fig.  7.  Simulated  degree  of  dissociation,  A*,  vs. 
equilibrium  temperature  compared  with  experimental 
fit  in  At  relaxation. 
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Timesups 

Fig.  8.  Percent  concentration  per  speciea  in  ther- 
mochemical  relaxation  of  air,  superheated  in  transla¬ 
tion  and  rotation  at  15000 if . 
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Normalized  Energy,  e,  and  Quantum  Number,  i 

Fig.  9.  Relax,  non  of  discrete  and  continuous  en¬ 
ergy  distribution?  toward  equilibrium  at  steady-state 
via  mixing  algorithm,  with  and  without  rotational- 
relaxation  steps. 
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